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SECTION  1 


INTRODUCTION 


1.1  OBJECTIVE 

The  main  objective  of  this  work  was  to  investigate  the  effects  of 
high  frequency  spiking  characteristics  of  the  High  Explosive  Simulation 
Technique  (HEST)  upon  testbed  and  test  article  response.  Spiking  is  the 
high  pressure,  high  frequency  deviation  of  a  HEST  loading  from  a  design 
nuclear  airblast  waveform.  This  investigation  provided  a  qualitative  and 
quantitative  evaluation  of  how  well  the  loading  from  a  HEST  cavity 
simulates  the  idealized  airblast  overpressure  of  a  nuclear  detonation.  A 
secondary  obj'ective  was  to  determine  how  the  frequency  content  of  a  HEST 
aata  record  influences  measured  structure/soil  response  date. 

1.2  BACKGROUND 

HEST  cavities  are  designed  to  match  Speicher- Brode  representations  of 
airblast  overpressure-time  waveforms,  but  significant  differences  often 
occur  between  the  recorded  data  and  predictions.  Due  to  the  discrete, 
rather  than  continuous  application  of  pressure  from  the  det-cord  and  the 
numerous  reflections  which  occur  within  a  HEST  cavity,  a  HEST  record  is 
strongly  characterized  by  high  frequency  (greater  than  1000  Hz)  and 
sometimes  high  magnitude  pressure  spikes.  A  HEST  record  can  also 
experience  lower  frequency,  low  magnitude  deviations  from  the  design 
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Spe’cher-Brode  pulsej  but  this  anomaly  is  less  common.  It  has  long  been 
of  concern  how  these  HEbi  deviations  from  an  ideal  airblast  time  history 
affect  soil  media  and  structure  response  within  a  testbed.  Is  the  high 
frequency  content  of  a  HEST  record  filtered  out  as  the  load  transfers 
through  soil  media  and/or  a  structure?  Or  does  the  high  frequency  input 
content  excite  high  frequency  response  modes?  These  questions  are 
difficult  to  answer  in  the  time  domain.  Frequency  domain  analysis  tools 
such  as  the  Fast  Fourier  Transform  (FFT),  Frequency  Response  Function 
(FRF),  and  the  inverse  FFT  are  ideal  for  investigating  these  questions, 
and  were  used  in  this  effort  to  obtain  a  qualitative  and  quantitative 
determination  of  HEST  fidelity  in  simulating  a  Speicher-Brode  nuclear 
environment. 


1.  Speicher,  S.J.  and  Brode,  H.L.,  Airblast  Overpressure  Analytic 

Expressions  for  Burst  Height,  Range  and  Time--0ver  an  Ideal  Surface, 
PSR  Note  38^,  Pacific-Sierra  Research  Corp.,  Los  Angeles,  CA, 
November  1981,  as  modified  for  time  of  arrival  at  high  overpressures 
by  memo  from  S.J.  Speicher,  Pacific-Sierra  Research  Corp.,  7  June 
1982. 


SECTION  2 


THEORETICAL  DEVELOPMENT 
2.1  CAUSE-EFFECT  RELATIONSHIP 

The  fundamental  assumption  in  the  analysis  v/hich  follows  is  that 
there  is  a  linear  cause-effect  relationship  between  two  sets  of 
corresponding  data  records.  A  linear  cause-effect  relationship 
necessarily  means  that  there  is  an  input  and  an  output,  related  by  a 
linear  transfer  function  (see  Figure  1).  The  input,  x(t),  may  be  a  HEST 
overpressure-time  waveform,  and  the  output,  a(t),  may  be  a  soil  or 
structure  response-time  waveform.  Any  test  data  waveform  pair,  x(t)  and 
a(t),  will  contain  a  certain  amount  of  data  acquisition  system  noise. 

This  analysis  will  ignore  the  presence  of  noise  and  assume  that  x(t)  and 
a(t)  are  actual  loading  and  response  behavior.  This  analysis  procedure 
should,  therefore,  not  utilize  records  containing  high  levels  of  noise. 
Steps  should  be  taken  to  remove  noise  content  or,  if  this  is  not  possible, 
the  record  should  be  discarded. 

The  actual  transfer  function  betv>een  an  input  data  record  and  an 
output  data  record  may  be  highly  nonlinear.  IJonlinearities  may  result 
from  concrete  crushing  and  cracking,  steel  yielding,  and  soil  media 
deforming  inelastical ly ,  or  from  undefined  soil -silo  interactions. 

However,  a  linear  assumption  may  be  acceptable  when  investigating  changes 
in  output  resulting  from  slight,  variations  in  input,  such  as  might  be  the 
case  when  considering  the  di-’^erence  between  a  HEST  record  and  a  best-fit 


3 


Speicher-Brode  nuclear  airblast  waveform.  If  there  are  significant 
variations  between  different  input  data  records,  the  assumed  linear 
transfer  function  may  cause  some  alterations  in  output  which  are 
erroneous.  The  amount  of  variation  in  the  input  which  can  be  allowed, 
without  causing  significant  erroneous  alterations  in  the  output,  is 
currently  indeterminate. 

2.2  FAST  FOURIER  TRANSFORM 

The  Fourier  transform  is  a  way  of  representing  a  time  domain 
function,  x(t),  in  the  frequency  domain.  If  X{(j)  is  the  Fourier  transform 
of  x{t),  then  x(t)  and  X((j)  are  called  a  Fourier  transform  pair. 

x(t)  <=>  X((o)  (la) 


X(<o)  =  n(t)e"''‘^^dt 

J  -<a 


(lb) 


where  t  =  time  (sec) 

Q  =  frequency  (rad/sec) 

1.VT 

If  x(t)  is  a  real  continuous  function  of  infinite  duration,  then  X(u)  is  a 
set  of  complex  numbers 

X(u)  =  a(u)  +  ib(«)  =  [A(«)V^^“'  (2a) 

which  define  both  the  amplitude  and  phase  associated  with  each  frequency, 
u,  of  the  function  x(t).  The  amplitude.  A,  and  phase  angle,  ?),  associated 
with  each  point  of  X(o)  are: 

A(u)  =  +  Cb(u)]2  =  [x(u)|  (2b) 

{5(u)  =  tan"^  (2c) 

ai  u) 

The  amplitudes.  A,  are  all  real  numbers  and  the  resulting  Fourier 
transform  amplitude  is  a  real  function  which  will  be  referred  to  as  |x(<j)j. 
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One  time  domain  function  of  particular  concern  is  the  dc  component,  which 
is  constant  with  time  (see  Figure  2).  The  Fourier  transform  of  a  dc 
component  is  a  single  value  at  zero  frequency,  known  as  a  Dirac-delta 
function.  The  significance  of  this  Fourier  transform  pair  will  be 
addressed  in  the  next  subsection. 

The  analysis  presented  in  this  report  utilizes  real  test  data 
records.  These  records  are  of  finite  duration  and  are  finely  digitized 
sets  of  points  with  a  constant  time  step.  Therefore,  this  analysis  must 
use  the  Discrete  Fourier  Transform  (DFT): 

X(n/NT)  =  4  x(KT)e‘‘'^’'"'^^^  (3) 

^  K=0 


where  x{KT) 
K 
n 


discrete  time  series 

0,  1,  2,  ...,  (N-1);  time  domain  counter 

0,  _^1,  +  2,  ...;  frequency  domain  counter 


T  =  time  step 

N  =  total  number  of  time  steps 
=  Tq  =  signal  duration 

The  DFT  has  to  be  scaled  to  approximate  the  continuous  integral  transform: 


X^(n/NT)  =  TpX(n/NT)  (4) 

where  X^(n/NT)  =  scaled  DFT. 

The  Fast  Fourier  Transform  (FFT)  is  a  computationally  efficient 
algorithm  for  calculating  the  DFT.  The  FFT  algorithm  reduces  the 
computation  of  the  DFT  of  an  N  point  time  series  from  operations  to 
Nlog^N  operations.  If  N  is  ';;rly  large  the  savings  in  the  number  of 
operations  can  be  very  significant.  For  example,  if  the  number  of  points 
in  a  time  series  is  4096,  Equation  (3)  requires  16,777,216  operations, 
while  the  FFT  requires  only  49,152. 


5 


The  number  of  points  in  the  FFT  array  is  one  half  the  number  of 
points  in  the  original  time  series.  The  first  FFT  value  is  defined  at 
zero  frequency;  the  second  FFT  value  is  defined  at  the  fundamental 
frequency: 

f  = -  1 - ^  (5) 

0  record  duration 

and  the  third  FFT  value  is  defined  at  twice  the  fundamental  frequency, 

2fQ,  and  so  on.  The  final  FFT  value  is  defined  at  the  Nyquist  frequency:. 

f  _ _ _ r-  (6) 

c  2 (data  record  time  step) 

The  Nyquist  frequency,  f^,  should  be  greater  than  or  equal  to  the 
highest  frequency  of  concern  in  any  subsequent  analysis. 

One  further  step  is  application  of  a  "weighting"  function  or  "window" 
to  a  discrete  timp  sei^ies  before  it  is  transformed.  It  is  best  that  a 
discrete  time  series  start  and  end  up  at  zero  in  order  to  prevent 
"leakage"  (Iruncation)  effects.  Leakage  is  attenuation  i.<  amplitude  of 
the  primary  frequency  component,  and  magnification  of  ether  secondary 
frequency  components.  "Weighting"  functions  or  "windows"  are  use1  to 
gradually  bring  the  beginning  and  end  portions  of  a  data  trace  back  to 
zero.  Most  test  data  traces  start  at  zero,  but  most  do  not  return  to  zero 
at  the  end.  For  the  analysis  in  this  report,  a  cosine  squared  spline  was 
applied  to  the  final  15  percent  of  all  data  traces  to  return  them  to  zero. 
2.3  INVERSE  FAST  FOURIER  TRANSFORM 

The  original  discrete  time  series,  x{KT),  can  be  reconstructed  from 
the  inverse  of  the  OFT. 

x(KT)  =  X(n/NT)e^'^’^"'^/^^  (7) 

n=0 
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Again,  an  inverse  FFT  algorithm  is  used  for  computational  efficiency.  The 
new  resultant  time  series  differs  slightly  from  the  original  time  series. 
As  an  example,  notice  the  three  plots  in  Figure  3.  An  original  discrete 
time  series  is  shown  in  Figure  3a.  The  FFT  amplitude  spectrum  is  shown  in 
Figure  3b  and  the  inverse  FFT  is  shown  in  Figure  3c.  If  there  v;ere  no 
inaccuracy  in  the  FFT  algorithm,  the  plots  in  Figures  3a  and  3c  would  be 
identical.  An  obvious  difference  is  that  the  inverse  FFT  appears  to  be 
offset  from  the  time  axis.  The  offset  appears  to  be  constant  with  time, 
since  the  final  portion  of  the  inverse  FFT  ends  up  at  the  same  level  as 
tne  beginning.  A  constant  offset  with  time  indicates  the  presence  of  a  dc 
component  in  the  inverse  FFT  time  history.  As  was  discussed  in  the 
previous  subsection  and  shown  in  Figure  2,  a  dc  component  in  a  time 
history  is  caused  by  the  value  of  the  FFT  at  zero  frequency.  The  FFT  in 
Figure  3b  does  in  fact  have  a  value  at  zero  frequency.  But  zero  frequency 
is  related  to  an  infinite  duration  in  the  time  domain  and  the  original 
discrete  time  series  in  Figure  3a  is  of  finite  duration.  Due  to  the 
finite  time  duration,  the  value  of  the  FFT  at  zero  frequency  must  be 
erroneous,  and  is  therefore  the  reason  for  the  offset  in  Figure  3c.  If 
the  dc  component  is  removed  from  the  inverse  FFT,  it  becomes  nearly 
identical  to  the  original  discrete  time  series  (see  Figure  4).  The 
inverse  FFT  removes  some  of  the  high  frequency  spiking  present  in  the 
original  time  history  due  to  the  reduced  number  of  points  defining  the  new 
record.  The  value  of  the  offset  is  equal  to  tv;enty  times  the  value  of  the 
full  integral  of  the  original  discrete  time  series.  The  absolute  value  of 
the  full  integral  of  the  original  discrete  time  series  is  also  equal  to 
the  value  of  the  FFT  amplitude  spectrum  at  zero  frequency. 
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2  A  FREQUENCY  RESPONSE  FUNCTION 


The  formuletion  of  ;;  frequency  response  function  (FRF)  requires  an 
input-output  re1 ationrnip  described  in  Section  2.1.  An  FRF  is  a  transfer 
function  in  the  frequency  domain  which  completely  defines  the  dynamic 
characteristics  of  a  lineav  system.  Again,  using  the  same  nomenclature  as 
was  used  in  Section  2.1,  x(t)  represents  an  input  time  history  and  a(t) 
represents  an  output  time  history.  The  Fourier  trarisfonns  of  x(t)  and 
a{t)  (X{o)  and  A(u),  respectively)  are  a  set  of  complex  numbers.  The  FRF, 
H(u),  is  simply  the  Fourier  transform  of  the  output  divided  by  the  Fourier 
transform  of  the  input,  and  is  also  a  set  of  complex  numbers. 

H(<o)  =  A(,o)/X(u)  (8) 

At  each  individual  frequency  an  FRF  describes  the  output  response  of  a 
linear  system  subjected  to  an  input  defined  by  a  crnstant  amplitude  sine 
wave  of  fixed  frequency.  The  input  is  of  the  form: 

x(t)  =  x^sin  yt  (9) 

The  output  response  will  be  a  sine  wave  at  the  same  frequency,  y,  fixed 
amplitude,  a^,  and  phase  difference,  /4; 

a(t)  =  aQsin(yt  -  i^)  (10) 

From  Reference  2: 

Information  about  the  amplitude  ratio  a  /x  and  the  phase  angle  i 
defines  the  transmission  characteristics  or  transfer  function  of  th-^ 
system  at  the  fixed  frequency  y.  The  FRF  H(y)  results  if  the  amplitude 
ratio  and  phase  angle  can  be  plotted  as  a  function  of  frequency... 

Instead  of  thinking  of  amplitude  ratio  and  phase  angle  as  two  separate 
quantities,  it  has  become  customary  in  vibration  theory  to  use  a  single 
complex  number  to- represent  both  quantities.  This  is  H(y)  which  is 
defined  so  that  its  magnitude  is  equal  to  the  amplitude  ratio  and  the 
ratio  of  its  imaginary  part  to  its  real  part  is  equal  to  the  tangent  of 
the  phase  angle. 


2.  Stearns,  S.D.,  Digital  Signal  Analysis,  Hayden  Book  Co., 
Rochelle  Park,  NJ,  1975. 
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If 


H(d)  =  B(u)  +  iC(u)  (11a) 

then 

|h(u)|  +  C2  =  ao(u)/xo  (lib) 

«5(u)  =  tan“^  ^  (11c) 

The  FRF  amplitude  ratio,  jH(u)j,  is  a  useful  tool  in  determining  how  much 
of  the  HEST  pressure  records  gets  through  to  soil  and  structural  response 
at  a  given  frequency. 

2.5  FOURFIT 

Reference  3  discusses  the  purpose,  use  and  theory  of  the  program 
FOURFIT.  It  is  basically  a  program  which  will  examine  a  HEST  pressure 
record  and  determine  the  yield  and  peak  overpressure  of  a  "best-fit" 
Speicher-Brode  ideal  nuclear  airblast  waveform,  using  frequency  domain 
analysis.  The  ideal  airblast  waveform  can  be  substituted  for  a  HEST 
pressure  record  as  input  in  order  to  determine  new  response  data.  A 
listing  of  FOURFIT  is  presented  in  Appendix  A  with  slight  modifications 
for  saving  the  "best-fit"  airblast  waveform. 

2.6  MODIFIED  OUTPUT  RESPONSE 

Since  the  FRF  is  a  linear  transfer  function  which  completely 
describes  the  dynamic  characteristics  of  a  linear  system,  varied  input 
waveforms  can  be  applied  to  the  system  to  get  new  outputs.  As  long  as 
changes  in  the  input  are  relatively  minor,  the  linear  assumption  remains 


3.  Steedman,  D.W.  and  Partch,  J.C.,  F0URFIT--A  Computer  Code  for 
Detennining  Equivalent  Nuclear  Yi¥ld  and  Peak  Overpressure  by  a 
Fourier  Spectrum  Fit  Method,  as  yet  unpublished  DNA  report,  Appl i ed 
Research  Associates,  Inc.,  Albuquerque,  NM,  May  1984. 
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valid  (see  Section  2.1).  Repeating  Equation  (8)  the  FRF  H((j)  is  defined 
as: 

H(u)  =  A((j)/X((j)  (8) 

Rearranging  this  equation  one  obtains: 

A{u)  =  H(u)*X(u)  (12) 

The  original  input  time  history  is  x(t).  A  new  input  time  history  x'(t) 
can  be  Fourier  transformed  to  obtain  X'(uj.  This  new  Fourier  transform 
can  be  substituted  into  Equation  (12)  to  obtain  a  new  output  Fourier 
transform,  A'(u): 

A'{o)  =  H(u)*X'(-o)  (13) 

A'((j)  can  be  inverse  Fourier  transformed  to  obtain  a  new  output  response 
time  history,  a' (t) . 
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SECTION  3 


PROGRAM  RESULTS 

3.1  PROGRAM  FREQRES 

Program  FREQRES  calculates  an  FRF  for  a  given  pair  of  input  and 
output  time  histories.  Using  the  defined  FRF,  the  program  will  calculate 
a  new  output  time  history  for  any  new  input  time  history  the  user  wants  to 
specify.  A  user's  manual  for  FREQRES  is  presented  in  Appendix  B,  and  a 
listing  of  the  program  is  presented  in  Appendix  C. 

3.2  TEST  DATA 

Data  from  a  HEST  test  of  a  surface  flush  vertical  silo  surrounded  by 
soil  were  analyzed  using  the  program  FREQRES.  Both  the  silo  and  the 
surrounding  soil  were  loaded  by  the  HEST  cavity.  The  test  data  had  a 
duration  of  49.85  msec  and  were  recorded  at  a  digitizing  rate  of 
200, OOQ  Hz.  This  resulted  in  a  time  step  of  5  x  10“^  sec,  and  a  total 
of  9970  digitized  points.  Figures  5  through  22  present  plots  of  the  18 
data  records  used  in  this  analysis.  Table  I  presents  a  brief  description 
of  each  of  the  plots. 

3.3  FOURFIT  RESULTS 

A  FOURFIT  analysis  was  performed  on  the  two  HEST  pressure  time 
histories,  test  data  record  number  2  on  the  structure  and  record  number  4 
on  the  free  field  soil.  Three  plots  result  from  a  FOURFIT  calculation: 

(1)  a  frequency  domain  plot  comparing  the  FFT  of  the  test  data  to  the  FFT 
of  the  "best-fit"  Speicher-Brode  ideal  nuclear  v/aveform,  (2)  a  time 
history  comparison  of  the  pressure  waveforms,  and  (3)  a  time  history 
comparison  of  the  resulting  impulse  curves.  Figures  23a,  23b,  and  23c 
present  these  three  plots  for  HEST  record  number  2  on  the  structure  and 
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Figures  24a,  24b,  and  24c  present  the  three  plots  for  HEST  record  number  4 
on  the  free  field  soil.  The  "best-fit"  Speicher-Brode  ideal  nuclear 
waveform  to  HEST  record  number  2  is  a  19.08  kt  yield  surface  burst  with  a 
peak  overpressure  of  20,280  psi.  For  HEST  record  number  4  it  is  a  7.83  kt 
yield  surface  burst  with  a  peak  overpressure  of  14,620  psi.  The  above 
results  indicate  substantial  variation  in  the  effective  yield  and  peak 
overpressure  across  the  testbed,  and  between  measurements  on  the  structure 
and  on  free  field  soil. 

3.4  FREQUENCY  RESPONSE  FUNCTIONS 

If  the  system  represented  by  the  FRF  transfer  function  is  linear,  the 
absolute  magnitude  of  the  FRF  amplitude  ratio  has  significance.  For 
example,  if  one  specifies  a  pressure  record  in  psi  as  an  input  data  record 
and  a  velocity  time  history  in  in/s  as  an  output  data  record,  then  the 
amplitude  ratios  of  the  FRF  can  be  multiplied  by  a  constant  pc  (p  is  the 
density  of  the  material  in  the  linear  system  and  c  is  the  loading  wave 
velocity)  to  normalize  the  FRF.  In  a  normalized  amplitude  ratio  FRF,  a 
value  of  1.0  at  a  particular  frequency  represents  perfect  transmission  of 
pov/er  at  that  frequency.  A  value  of  less  than  1.0  represents  a  decay  in 
power  and  a  value  greater  than  1.0  represents  an  amplification  of  power. 

In  a  reinforced  concrete  silo  pc  is  fairly  constant,  and  therefore  can  be 
used  to  normalize  the  FRF.  For  soil  response  pC  may  vary  significantly 
with  depth  and  time  and,  therefore, a  constant  does  not  exist  to  normalize 
the  FRF.  For  other  input/output  combinations,  different  constants  of 
proportionality  exists.  If  the  input  is  pressure  and  the  output  is 
strain,  then  the  constant  will  be  a  stiffness  modulus.  If  the  input  is 
pressure  and  the  output  is  stress,  then  no  constant  of  proportionality  is 
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necessary  since  the  input  and  output  are  already  in  the  same  units  and  the 
FRF  is  already  normalized. 

Figure  25  shows  two  FRF's  for  vertical  soil  stress  response  at  0.5  ft 
and  5.21  ft  depths.  Notice  that  at  the  0.5'  depth,  power  transmission 
from  the  HEST  pressure  loading  to  the  soil  stress  response  is  strong  at 
most  frequencies  from  0  to  3000  Hz.  In  fact  there  is  significant  power 
transmission  at  1700  Hz,  2200  Hz,  and  2800  Hz  (relatively  high 
frequencies).  The  reason  for  this  could  be:  (1)  noise  in  the  data  record 
at  the  three  frequencies  mentioned  above,  or  (2)  exitation  of  natural 
frequencies  in  response.  Notice  from  Figure  23a  that  the  input  power 
content  at  frequencies  betv/een  1700  Hz  and  2800  Hz  is  very  low  compared  to 
input  power  at  frequencies  less  than  100  Hz.  Since  the  HEST  contains  low 
power  at  the  higher  frequencies,  strong  power  transmission  at  these 
frequencies  still  results  in  relatively  low  power  for  the  FFT  amplitude 
spectrum  of  soil  stress  response  (see  Figure  3).  Also  note  from  Figure  25 
that  power  transmission  at  all  frequencies  greater  than  1000  Hz  is 
dramatically  reduced  when  going  from  the  0.5  ft  depth  to  the  5.21  ft 
depth.  The  power  transmission  between  100  and  1000  Hz  is  also 
significantly  reduced,  but  not  to  the  degree  evident  at  higher 
frequencies.  This  suggests  that  by  5  ft  depth,  the  soil  has  significantly 
filtered  out  the  high  frequency  characteristics  of  the  HEST  pressure 
loadi ng. 

Figure  26  shov/s  two  FRF's  for  vertical  soil  velocity  response  at 
5.21  ft  and  12.21  ft  depths.  The  FRF  at  the  5.21  ft  depth  is  very  similar 
in  shape  to  the  FRF  at  the  same  depth  in  Figure  25.  There  is  a  constant 
decay  in  power  transmission  between  100  Hz  and  1000  Hz.  At  the 
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12.21  ft  depth,  power  transmission  is  even  further  reduced  between  100  and 
1000  Hz.  At  this  depth  power  transmission  is  relatively  low  at  all 
frequencies  above  350  Hz.  Soil  is  a  very  good  attenuator  of  the  high 
frequency  power  of  a  HEST  cavity,  starting  with  the  higher  frequencies. 

Figure  27  shows  two  FRF's  for  vertical  structural  velocity  response 
at  0.83  ft  and  3.28  ft  depths.  The  FRF  amplitude  ratio  can  be  normalized 
through  application  of  the  proportionality  constant,  pC.  Assume  the 
density  cf  concrete  to  be  4.7  slugs/ft  and  the  loading  wave  velocity  to 
be  10,000  ft/s. 

pC  =  (4.7  si ugs/ft^)(  10,000  ft/s){f rVl728  in^) 

,  (14) 

=  27.2  Ib-s/in*^ 

Multiplying  this  constant  times  the  amplitude  ratio  scale  of  0.0  to  0.03 
in  Figure  27  results  in  a  normalized  scale  of  0.0  to  0,82.  A  noticeable 
large  frequency  powe*^  transmission  (75  percent  on  the  normalized  scale) 
exists  at  approximately  250  Hz.  This  is  related  to  the  natural  frequency 
of  axial  response  of  the  vertical  cylinder  test  article.  The  structure 
was  approximately  21.5  ft  in  length.  At  the  0.83  ft  depth  the  time  it 
took  the  axial  stress  wave  to  reach  the  bottom  of  the  cylinder  and  reflect 
back  up  is  (assuming  shock  wave  velocity  in  concrete  =  10,000  ft/s): 


2(21.5  ft  -  0.83  ft) 
1U,UUU  tt/s 


4.13  x  10"^  sec 


(15) 


The  frequency  associated  with  4.13  x  10"*  sec  is  1/4.13  x  10"  sec 
which  equals  242  Hz.  At  the  3.28  depth  the  travel  time  and  associated 
frequency  ar’  3.64  x  10"^  sec  and  274  Hz,  respectively.  The  FRF's  in 
Figure  27  also  show  a  decay  in  power  transmission  with  depth,  similar  to 
that  for  soil  response.  Comparing  Figures  25  and  27,  the  decay  with  depth 
is  not  as  dramatic  in  the  structure  as  it  is  in  the  soil. 


The  FRF's  for  structural  axial  strains  and  hoop  strains  (see  Figures 
28  and  29)  shov/  no  clear  pattern  of  power  transmission  decay  with  depth. 
This  indicates  that  the  structure  tends  to  transmit  most  of  the  power  from 
a  REST  pressure  time  waveform  over  the  broad  frequency  range  of  0  to 
3000  Hz.  The  normalizing  proportionality  constant  for  structural  strains 
is  assumed  to  be  the  constrained  modulus  for  concrete*  since  the  concrete 
in  the  structure  is  confined  with  a  high  percentage  of  steel.  Assume  a 
constrained  modulus  of  5.16  x  10°psi.  The  strain  test  data  are  in 
units  of  micro-strain,  such  that  1  x  10^  has  to  be  factored  out  of  the 
proportionality  constant.  The  resulting  constant  is  5.16.  Multiplying 
this  constant  times  the  amplitude  ratio  scales  of  0.0  to  0.8  in  Figures  28 
and  29  results  in  a  normalized  scale  of  0.0  to  4.13.  A  normalized 
amplitude  ratio  of  1.0  corresponds  approximately  to  0.2  on  the  scales  in 
Figures  28  and  29.  There  are  two  characteristics  common  to  both  Figures 
28  and  29.  At  the  5.2  ft/5.3  ft  depth  there  is  strong  pov/er  transmission 
(335  per  cent  on  the  normalized  scale)  at  150-200  Hz.  At  the  1.3  ft  depth 
there  is  also  strong  pov/er  transmission  (181  percent  on  the  normalized 
scale)  at  1200-1600  Hz.  Why  the  150-200  Hz  strong  power  transmission  is 
peculiar  to  the  o.Z  ft/5.3  ft  depth  is  currently  unclear.  The  only 
structural  response  mode  with  a  natural  frequency  as  low  as  150-200  Hz  is 
the  axial  response  of  the  cylinder  associated  with  the  axial  stress  wave 
traveling  back  and  forth  down  the  entire  length  of  the  cylinder  and 
reflecting  off  each  end.  But  if  this  were  the  cause  of  the  strong 
frequency  power  transmission  at  the  5.2  ft/5.3  ft  depth,  then  it  should 
also  occur  at  the  other  depths  as  well.  It  does  not.  The  1.3  ft  depth  in 
the  cylinder  occurs  in  a  thick  walled  portion  of  the  cylinder  called  the 
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headworks.  Hoop  expansion  of  the  cylinder  due  to  passage  of  the  axial 
compressive  wave,  and  also  hoop  compression  due  to  large  ground  shock 
stresses  surrounding  the  cylinder  are  very  strong  near  the  surface.  Hoop 
expansion  and  compression  are  associated  with  the  breai.hing  mode  response 
of  a  cylinder.  The  natural  frequencies  for  the  breathing  mode  response  of 
a  cylinder  are  (Ref.  4,  pg.  298,  Table  12-1): 


f. 


1 


(16) 


where  f^.  =  ith  natural  frequency 

i  =  response  mode  (=0  for  breathing  mode) 

E  =  Young's  Modulus 
R  =  cylinder  radius  to  midsurface 
It  =  density  of  shell  material 
V  =  Poisson's  ratio 
=  (1  +  i^)^^^  =  1  for  i  =  0 

The  following  values  are  substituted  into  Equation  (16)  to  determine  the 
lowest  breathing  mode  natural  frequency  of  the  headworks: 
i  =  0 

E  =  5.77  X  10^  psi  (effective  Young's  Modulus  including 
concrete  and  steel  contributions) 

R  =  20.0  inches 
n  =  2.25  X  10"^  lb  s^/in^ 

V  =  0.2 

Xo  =  1 

The  result  is  1300  Hz.  Notice  that  this  falls  within  the  1200-1600  Hz 
range  of  strong  power  'transmission  for  the  1.3  ft  depth  in  both  Figures  28 
and  29. 


4.  Blevins,  R.D.,  Formulas  for  Natural  Frequency  and  Mode  Shape, 
Van  Nostrand  Reinhold  Co.,  New  York,  NY,  1979. 
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3.5  FAST  FOURIER  TRANSFORM  COMPARISONS 


Figure  30  shows  a  comparison  of  the  original  FFT  amplitude  spectrum 
for  test  data  record  number  6  and  the  FFT  amplitude  spectrum  after  the 
Spei cher-Brode  waveform  influence  has  been  included  (see  Section  2-6). 

The  FFT  with  Speicher-Brode  influence  has  higher  power  at  most  frequencies 
between  100  and  3000  Hz.  Notice  that  the  same  is  true  when  comparing  the 
Speicher-Brode  "best-fit"  waveform  FFT  to  HEST  record  number  4  FFT  in 
Figure  24a.  In  fact  the  higher  power  in  the  Speicher-Brode  "best- fit" 
waveform  FFT  is  the  cause  of  the  higher  power  in  the  FFT  with 
Speicher-Brode  influence  in  Figure  30.  But  this  higher  power  has  very 
little  effect  upon  the  inverse  FFT  time  history,  as  will  be  illustrated  in 
the  next  section.  For  test  record  number  18  which  Uf^ed  HEST  record 
number  2  instead  of  ^  as  the  input  data  »"ecord  in  FREORES,  the  increased 
power  is  only  evident  above  800  Hz  (see  Figure  31).  Notice  that  the  same 
is  true  in  Figure  23a  in  which  the  fit  between  100  and  700  Hz  is  very  good. 
3.6  TIME  HISTORY  COMPARISONS 


For  records  5  through  20,  there  is  no  visible  difference  between  the 
original  record  time  histories  and  the  response  time  histories  had  the 
surface  pressure  loading  been  an  ideal  Speicher-Brode  nuclear  waveform. 

As  an  example.  Figure  32  presents  record  number  5  and  its  Speicher-Brode 
input  comparison  waveform.  The  high  frequency  spiking  characteristic  of 
HEST  pressure  records  has  negligible  effect  upon  soil  and  structure 
response. 
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SECTION  4 


CONCLUSIONS  AND  RECOMMENDATIONS 

Analysis  of  the  test  data  presented  in  this  report  indicates  that  the 
high  frequency  spiking  characteristic  of  HEST  pressure  records  has 
negligible  effect  upon  testbed  and  test  article  response  when  compared  to 
loading  from  a  "best-fit"  ideal  Speicher-Brode  nuclear  airblast  v/aveform. 
The  high  frequency  power  content  of  a  HEST  pressure  loading  attenuates 
with  depth  in  both  soil  and  structure.  The  power  decay  is  more  dramatic 
in  soil  and  is  evident  in  both  stress  and  motion  response.  Structural 
strain  records  show  little  power  transmission  decay  with  depth.  This 
indicates  that  structural  strain  response  tends  to  transmit  most  of  the 
power  content  of  HEST  pressure-time  waveforms.  Beyond  the  2000-3000  Hz 
range  the  power  content  in  a  HEST  pressure-time  waveform  is  very  low 
compared  to  that  in  the  lower  frequencies,  so  that  power  transmission  at 
the  higher  frequencies  is  insignificant. 

Strong  power  transmission  is  evident  in  structural  response  FRF's  at 
natural  response  mode  frequencies.  Natural  modes  are  excited  by  the  broad 
frequency  range  of  power  content  of  a  HEST  pressure  loading,  such  that  the 
FRF  amplitude  ratio  becomes  magnified  at  the  modal  frequencies. 

The  HEST  pressure  records  used  in  this  report  (Figures  5  and  6)  were 
very  good  pressure  waveforms,  in  which  the  Speicher-Brode  "best-fits"  from 
the  FOURFIT  program  matched  the  pressure  and  impulse-time  waveforms  very 
closely  (see  Figures  23b,  23c,  24b,  24c).  It  would  be  interesting  to  run 
this  analysis  with  HEST  records  from  a  test  which  do  not  provide  such  good 
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representations  of  ideal  Speicher-Brode  nuclear  airblast  v^aveforms.  In 
this  case  there  might  be  greater  variation  in  output  response 
comparisons.  Also,  another  interesting  use  of  this  analysis  procedure 
would  be  to  determine  what  testbed  and  test  article  response  would  have 
been, had  the  HEST  loading  been  exactly  as  original  design.  Instead  of 
using  the  "best-fit"  Speicher-Brode  waveform  to  a  HEST  record,  have  the 
FREQRES  program  read  in  the  original  design  waveform. 

The  analysis  procedure  outlined  in  this  report  can  be  used  to 
investigate  the  fidelity  of  any  nuclear  airblast  simulation  technique.  It 
is  not  restricted  to  HEST. 
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APPENDIX  A 


LISTING  OF  PROGRAM  FOURFIT 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


PROBRAM  FOURFIT (INPUT,  OUTPUT, TAPE5=INPUT,  TftPE6=0UTPUT, 
♦  TOPES,  TAPES6,  T0PE48,rT0PE49j) 

*-»*-)HHt-**-***********'JHHH.**»*********%'«"»****-)HHHH(-***** 

PROGROM  FOURFIT  ESTIMOTES  THE  PEAK  OVERPRESSURE 
AND  NUCLEAR  YIELD  FOR  AIRBLAST  SIMULATION  RECORDS 
BY  COMPARING  FITS  OF  THE  DATA  FOURIER  AMPLITUDE 
SPECTRUM  TO  THE  FOURIER  AMPLITUDE  SPECTRA  OF  TRIAL 
SPEICHER-BRODES.  RESULTS  ARE  WRITTEN  TO  A  FILE 
(TAPE4a)  TO  BE  READ  AND  PLOTTED  BY  PROGRAM  FOURPLT. 
***************4*********************************** 


COMMON' /FFT  /  FRQ (3001 ), AMP (3001 ),  XFFT (3001 ) 

COMMON  /ITERAT/  W (5) , P (5) , DELTAW (5) , DELTAP (5) , YLD (5) 

COMMON  /THIST  /  TTIM (6000) , PRESS ( 18000) , TIMP (S999) , PIMP (S9S9) , 
*  PFILT(6000) 

COMMON  /IMP  /  IIMP,DTD, DTB, TPEB, DTBN  ■ 

COMMON  /POINTS/  NEPTS, NBPTS, NI, NEF, NBF 
COMMON  /ESTIM  /  PSOI , WI , PP, W13, PSOF, WF, FS0 
COMMON  /PEAK  /  DP, TP, PSO, ALPF 
COMMON  /SBCONS/  RSKFT, YS, S, XM 

COMMON  /FILT  /  IFILT, FLO (7) , PFDMX (7) , PFBMX (7) 

COMMON  /PLOTV  /  ITL (8) , ISTL (8) , IDB 
COMMON  /UNITS  /  lUNITS, JUNITS 
COMMON  /COUNT  /  ICOUNT, lOPT, LFILT 
COMPLEX  XFFT 


TAPES  CONTAINS  INPUT  PARAMETERS 

NEPTS=  NO.  OF  POINTS  TO  BE  READ  FROM  TAPE 
IUNITS=1  FOR  TAPE  INPUT  PRESSURE  IN  PSI 
=-l  FOR  TAPE  INPUT  PRESSURE  IN  MPA 
JUNITS=1  FOR  tpipe  INPUT  TIME  IN  MILLISECONDS 
=-l  FOR  TAPE  INPUT  IN  SECONDS 
PSOI=INITIAL  PEAK  OVERPRESSURE  ESTIMATE  IN  MPA 
WI=INITIAL  NUCLEAR  YIELD  ESTIMATE  IN  KT 
IOPT=l;  FITTING  ROUTINE  TO  BE  DONE 
IOPT=S:  JUST  FOURIER  TRANSFORM  THE  DATA 
IOPT=3:  JUST  FOURIER  TRANSFORM  THE  SPEICHER-BRODE 
DEFINED  BY  PSOI,  WI 
IFILT=1  FOR  FILTER  TO  BE  EXECUTED 
IFILT=-1  FOR  NO  FILTER 

FLO=LOW  END  CUTOFF  FREQUENCY  (UP  TO  7  ALLOWED) 
(NOTE; FOR  LESS  THAN  7  FILTERS,  FLO 
MUST  BE  SET  TO  0.  TO  ESCAPE  THE  LOOP. ) 
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REWIND  £ 

REPID  (£,  111)  NEPTS,  lUNITS,  JUNITS 
REfiD<£, 11£)  PSOI,WI 
REPlD<£,  113)  IQPT,  IFILT 
REPlD<£,  115)  (FLOU),  1  =  1,7) 

111  FORMAT (315) 

11£  FORMAT (£F5. £) 

113  FORMAT (£15) 

115  FORMAT (7F10.0) 

WRITE(6, 1)  PSOI, WI 

1  FORMAT (£X,  *PSOI=*,  F5.  £,  5X,  *WI=*,  F5. £) 

WRITE(4a, 113)  lOPT, IFILT 

ICOUNT=0 

NBPTS=£0Aa 

IFdOPT.  EQ.  3)  eO  TO  7 

3  CALL  EBREAD 

IFdOPT.  EQ.  £)  GO  TO  666 

4  CALL  FIT 
7  ICOUNT=l 

CALL  RANGE 
CALL  SPBRODE 
666  END 

SUBROUTINE  EBREAD 

C  ********»***»********************»*********** 

C  THIS  SUBROUTINE  READS  PRESSURE  VALUES  FROM  AN 

C  EBCDIC  TAPE  BASED  UPON  THE  FORMAT  PREVIOUSLY 

C  USED  BY  WES. 

C  *****************»******«*********-»********** 

c 

c 

COMMON  /FFT  /  FRQ (3001 ), AMP (3001 ), XFFT (3001 ) 

COMMON  /POINTS/  NEPTS,  NBPTS,  NI ,  NEF, NBF 

COMMON  /THIST  /  TTIM (6000) ,  PRESS ( 1£000) ,  TIMP (£9S9) , PIMP (£999) , 
*  c.FILT(6000) 

COMMON  /FILT  /  IFILT, FLO (7) , PFDMX (7) , PFBMX (7) 

COMMON  /IMP  /  IIMP,DTD,  DTB,  TPE&,DTBN 
COMMON  /UNITS  /  lUNITS, JUNITS 
COMMON  /PLOTV  /  ITL (S) , ISTL (8) , IDB 
COMMON  /COUNT  /  ICOUNT, lOPT, LFILT 
COMPLEX  XFFT 

DIMENSION  IWKE (5500) ,  WKE (5500) 

EQUIVALENCE  (I WKE ( 1 ) , WKE ( 1 ) ) 

DELP  IS  THE  DATA  BASELINE  SHIFT.  BE 
SURE  THAT  IT  IS  IN  THE  PROPER  UNITS. 

DIMENSION  DUM(3),DA(5) 

DELP=0. 0 
REWIND£b 

READ  TAPE  HEADER  INFORMATION 


22 


n  n  n  non 


READ (£6, 30) 

ITL(3) , 

ITL (4) , 

* 

DUM ( 1 ) , 

DUM(£) , 

* 

ITL(l), 

ITL(£), 

DTD, NP 

30  FORMAT (3 (£010) , E15. a, 15) 

DTD=£. *DTD 
NP=£*NEPTS 

ITL(5)=10H  PRESSURE 
ITL (a)=10HHISTORY 
ITL(7)=10H 
ITL (3) =10H 

WRITE(Ha,35)  (ITL  <L),L=1,8) 

35  FORMAT (8010) 

DO  £0  I=1,NEPTS 
TTIM(I)=0. 

PRESS (I) =0. 

£0  CONTINUE 

IF(EQF(£&))  900,901 

SET  UP  DATA  UNITS  CONVERSIONS; 

MSEC  TO  SEC  AND  PSI  TO  MPA. 

901  IF(JUNITS. GE, 1)  DTD=DTD*. 001 
PFACT=. 006894757 
IFdUNITS.  LT.0)  PFACT=1. 

IP=1 

TIME=0. 

NLINE=-NP/5 
RLINE=FLCAT (NP) /5. 

IF (RLINE. GT. NLTNE)  NLINE=NLINE+1 

READ  PRESSURE  VALUES 

DO  40  J=1,NLINE 

READ (£6, 50)  (DA ( JJ) , JJ=1 , 5) 

50  FORMAT (5E16. 8) 

IF (EOF (£6))  900, 90£ 

90£  DO  60  K=l,5 
P=D0(K) 

PRESS ( IP) = (P*PFACT) -DELP 
IP=IP+1 

60  CONTINUE 

40  CONTINUE 
IM=0 
TIME=0. 

DO  11  M=£,NP,£ 

IM=IM+1 

TTIM(IM)=TIME 

PRESS ( IM) = (PRESS (M) +PRESS (M-1 ) ) /£. 
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TIME=TIME+DTD 
11  CONTINUE 

SPLINE  THE  END  OF  THE  DPTA  TO  ZERO  IN 
CASE  OF  A  TRUNCATED  RECORD 
TLAST=TTIM<NEPTS) 

CALL  SPLINE (TLAST,  NEPTS,  TTIM,  PRESS) 
PMAX=0. 


IF  IDPT=1,  FIND  THE  TIME  TO  DATA 
PEAK  TO  AID  IN  PHASING  THE  OVERLAYS. 

AID  IN  PHASING  OVERLAYS 
PMAX=0. 

DO  78  IK=1,NEPTS 

PmAX=AMAX1 (PMAX,  press < IK) ) 

IF(PMAX. EQ. PRESS (IK) )  TPEB=TTIM ( IK) 
78  CONTINUE 


REMOVE  BASELINE  CORRECTION  FOR  POINTS 
BEFORE  THE  ARRIVAL  OF  THE  SHOCK 
DO  77  M=l, NEPTS 

IF<TTIM(M).  GT.  TPEB)  GO  TO  990 
PRESS  <M) =PRESS  <M) +DELP 
77  CONTINUE 

GO  TO  990 
900  WRITE(6, 70) 

70  FORMAT ( 10X,  *END-OF-FILE  REACHED  EARLY*,///) 
990  CONTINUE 

IFdOPT.  NE.  £)  GO  TO  45 

CALL  FMAX (PRESS,  NEPTS,  YPMN, Y^MX ) 

CALL  FMAX (TTIM,  NEPTS, XPMN, XPMX) 

WRITE (48,  100)  NEPTS,  XPMN,  XPMX, YPMN, YPMX 
IFdFILT.LT.  0)  GO  TO  700 

CALL  FOR  FILTERS  TO  BE  EXECUTED 
CALL  FLOOP (TTIM,  PRESS,  DTD, NEPTS, PFILT) 
RETURN 

700  WRITE(48,  105)  (TTIM (K) ,  K=1 ,  NEPTS) 

WRITE (48, 105)  (PRESS (KL) , KL=1, NEPTS) 

100  FORMAT (15, 4E15. 8) 

105  FORMAT (10E15. 8) 

45  IIMP=1 


IMPULSE 

CALL  IMPULSE ( I  IMP,  DTD,  NEPTS, NI ) 
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IFdOPT.  NE.  £>  GO  TO  110 
ITL(5)=10H  IMPULSE  H 
ITL<6)=10HISTQRY 
WRITE<48, 115)  ITL(5> , ITL(6) 

CALL  FMPX (TIMP, NI, XIMN, XIMX) 

CALL  FMAX<PIMP, NI,YIMN, YIMX) 

WRITE(48, 100)  NI, XIMN, XIMX, YIMN, YIMX 
WRITE<48, 105)  <TIMP<IH) , 1H=1, NI) 

WRITE<48, 105)  <PIMP(JH), JH=1, NI) 

115  FORMAT <SA10) 

C 

C  FIND  THE  FOURIER  TRANSFORM  AND  CALCULATE  AMPLITUDE. 

C 

110  TTOT=DTD*NEPTS 
C  FREQUENCY  INCREMENT 

DFE=1. /TTOT 
FQE=0. 

C  FOURIER  TRANSFORM 

CALL  FFTRC (PRESS, NEPTS, XFFT, IWKE, WKE) 

XRE=REAL<XFFT(1) ) / (S*NEPTS) 

XIE=AIMAG<XFFT<1) ) / <£*NEPTS) 

FQE=FQE+DFE 

FRQ(1)=:FQE 

C  AMPLITUDE  SPECTRUM 

AMP(1)=SQRT(£. *(XRE*XRE+XIE*XIE) )*TTOT 

NEF=NEPTS/£+l 

DO  80  JK=£, NEF 

FQE=FQE+DFE 

FRQ(JK)=FQE 

XRE=REAL(XFFT(JK) ) /NEPTS 
XIE=AIMAG (XFFT ( JK) ) /NEPTS 
AMP(JK)=SQRT(XRE*XRE+XIE*XIE)*TTOT 
80  CONTINUE 
C 

IFdOPT.  NE.  £)  RETURN 

ITL(5)=10H  FOURIER  A 

ITL (6)=10HMPLITUDE  S 

ITL(7)=10HPECTRUM 

CALL  FMAX (FRQ, NEF, XFMN, XFMX) 

CALL  FMAX (AMP, NEF, YFMN, YFMX) 

WRITE (48, 117)  ITL(5) , ITL (6) , ITL(7) 

117  FORMAT (3A10) 

WRITE (48, 100)  NEF, XFMN,  XFMX, YFMN, YFMX 
WRITE(48,  105)  (FRQ(LI),LI  =  1,NEF) 

WRITE(48, 105)  (AMP(JI), JI=1,NEF) 

RETURN 

END 

SUBROUTINE  FIT 

C  ***•»•)(•■»**********■»***♦•»*»•*■»***■»•******■)(■******■»***■»****■)(•**•)(■** 

C  THIS  SUBROUTINE  ITERATES  ON  YIELD  WITHIN  ITERATIONS  ON 
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C  PE«K  PRESSURE.  ITS  ftIM  IS  TO  REDUCE  THE  SUM  OF  THE  SQUARES 

C  OF  THE  DIFFERENCE  BETWEEN  THE  DATA  AMPLITUDE  AT  F(I)  AND 

C  THE  ESTIMATED  SPEICHER-BRODE  AMPLITUDE  AT  F(I)  DIVIDED 

C  BY  F(I)  BASED  UPON  A  TOLERANCE  ON  PEAK  PRESSURE  AND  YIELD. 

C  END  RESULT  IS  A  FINAL  ESTIMATE  OF  PEAK  OVERPRESSURE 

C  <PSOF)  AND  YIELD  <WF).  ALSO,  AN  ESTIMATE  OF  THE  SDODNESS 

C  OF  FIT  <DELL)  IS  DETERMINED.  PRESSURE  IS  IN  MPA, 

C  YIELD  IS  IN  KT. 

C  *•»*****■»■»****■)(•*■*•#***#***■«■*******•)(■■)(•**•»*********#*■»(■**■»*■»**** 

C 

c 

COMMON  /POINTS/  NEPTS,  NBPTS,  NI, MEF, NBF 

COMMON  /ESTIM  /  PSOI , WI , PP, W13, PSOF, WF, FS0 

COMMON  /ITERAT/  W (5) ,  P (5) ,  DELTAW (5) , DELTAP (5) , YLD (5) 

COMMON  /FFT  /  FRQ <3001 ),  AMP (3001 ), XFFT (3001 ) 

COMMON  /PEAK  /  DP,  TA,  PSO,  ALPF 
DATA  TOL/.01/ 

C 

P(l)=. l*PSOI 
P(e)=. 4*PS0I 
P(3)=1.0*PSOI 
P(4)=4.-»PS0I 
P(S)=10. *PSOI 
JPRESS=0 

LOOP  ON  PRESSURE  TOLERANCE 

DO  100  JJ=1,50 

JPRESS=JPRESS+1 
JMIN=£ 

JMAX=4 

IF<JPRESS. NE. 1)  GO  TO  105 
JMIN=1 
JMAX=5 

LOOP  ON  PRESSURE 

105  DO  £00  II=JMIN, JMAX 
PP=P (II) 

JYLD=0 
W(1)=0. 1*WI 
W(£)=0.  4-«-WI 
W(3)=l. 0*WI 
W(4)=4. 0*WI 
W(5)=10.  *WI 

LOOP  ON  YIELD  TOLERNACE 

DO  £50  KK=1,50 

JYLD=JYLD+1 
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IMIN=£ 

IMPiX=A 

IF(JYLD.NE. 1)  GO  TO  £55 

IMIN=1 

IMftX=5 

LOOP  ON  YIELD 

£55  DO  300  LL=IMiN, IMAX 

W13=W(LL)**. 33333 
IF (LL. NE. 1)  GO  TO  £56 
CALL  RANGE 

DETERMINATION  OF  RESIDUALS 

£56  DELTAW(LL>=0. 

DO  350  LK=1,NEF 

FSCL=FRQ(LK)*W13 
IF<FRQ(LK).GT.7000. )  GO  TO  300 
IF<FSCL. LT. FS0)  GO  TO  350 
CALL  AMPALG<FSCL, BAMP) 
ftMPN=ALOG10  < AMP  <LK) ) 

BAMPN=ALQG 1 0 ( BAMP ) 

DF£=FRG!  (LK)  #FRQ  <LK) 

DELTAA=<AMPN-BAMPN) /FRQ(LK> 

I F  <  FRQ  <  LK  > .  LT.  1 000.  >  DELTAA=£.  ♦DELTAA 
IF(FRQ(LK).GT.5000.  .AND.  FRQ  (LK) .  LT.  7000.  ) 

*  DELTAA=£.  *DELTAA 

DELTAA=DELTAA*DELTAA 
DELTAW (LL) =DELTAW  <LL) +DELTAA 
350  CONTINUE 
300  CONTINUE 

RESET  YIELDS 

EPSW=ABS(W(5)-W(1 ) )*£.  /(W(5)+W(l ) ) 

IF(EPSW. LT. TOL)  GO  TO  360 
CALL  RESETW 
£50  CONTINUE 

WRITE (6, 1£50) 

1£50  FORMAT (£X,*FAILED  TO  CONVERGE  ON  YIELD*) 

STOP  14 
360  CONTINUE 

DWMIN=AMIN1 (DELTAW ( 1 ) ,  DELTAW (£) ,  DELTAW (3) ,  DELTAW (4) , DELTAW (5) ) 
DO  365  MM=1,5 

IF(DELTAW(MM).ED.  DWMIN)  KW=MM 
365  CONTINUE 

YLD(II)=W(KW) 

DELTAP (II) =DELTAW (KW) 

£00  CONTINUE 
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RESET  PRESSURES 

EPSP=0BS(P(5)-P(1) )*£.  /(P(5)+P(l>  > 

IF(EPSP. LT. TOL)  GO  TO  400 
CALL  RESETP 
100  CONTINUE 

WRITE<G, 1100) 

1100  F0RMAT(2X,  i^FAILED  TO  CONVERGE  ON  PEAK  PRESSURE*) 

STOP  10 

400  DPMIN=AMIN1 (DELTAP ( 1 ) ,  DELTAP (2) ,  DELTAP (3) ,  DELTAP (4) ,  DELTAP (5) ) 
DO  405  NN=1,5 

IF (DELTAP (NN) . EQ, DPMIN)  KP=NN 
405  CONTINUE 

W13=YLD(KP)**. 33333 
PP=P(KP) 

DELL=DELTAP<KC')  /NEF 

RETURN 

END 

SUBROUTINE  AMPALG (FSCL,  SAMP) 

************■»****•»*•»***#***■*■*■«■***•»♦****»******•»#*■)(•■)(■**»«■**■)(■** 
THIS  SUBROUTINE  ESTIMATES  THE  FOURIER  AM^'LITUDE  OF  THE  TRIAL 
PEAK  PRESSURE  AND  YIELD  BASED  UPON  A  FIT  TO  THE  SUITE 
OF  NORMALIZED  SPEICHER-BRODE  FOURIER  AMPLITUDE  SPECTRA. 

THE  ALGORITHM  USES  SCALED  FREQUENCY  OF  INTEREST  (FSCL), 

SCALED  FUNDAMENTAL  FREQUENCY  OF  THE  S-B  OF  CONCERN  (FS0) 

AND  THE  PEAK  OVERPRESSURE  (PP)  TO  CALCULATE  THE  SCALED 
AMPLITUDE.  THE  ALGORITHM  USES  PRESSURE  IN  MPA  AND  YIELD 
IN  KT.  THE  EQUATIONS  ARE  FOR  A  SURFACE  BURST  ONLY.  THEY  ARE 
VALID  FOR  ANY  YIELD  AND  FOR  PEAK  OVERPRESSURE  UP  TO  100MPA 


COMMON  /ESTIM/  PSOI , WI , PP, W13, PSOF, WF, FS0 

Al=. 17a8*PP** (-. 7£) * (FSCL** (-1.  *PP** (-.  103) ) ) 

A£=. 01474*PP** (-. 15) *  <FSCL/FS0) **  <~1. 75) 

A3=.  001 1*PP** (PP** (-.  £34) ) *(FSCL/FS0) **(-£.  15) 

A4=. 0013£*FSCL**(-. 547) 

A5=. 01034*PP** (-. 1 13) * ( 1. /FSCL) * (FSCL/FS0) **(-1.5) 
A6=. 00001 1*PP*». 77* (FSCL/FS0) ** (-7.  5) 

A7=.  0000S66*PF'**.  3*  ( FSCL/FS0)  **  (  —  1 . 5) 

ASCL= A 1 - A£+ A3+ A4+A5-A6+A7 

BAMP=ASCL*PP*W13 

RETURN 

END 

SUBROUTINE  RESETW 

THIS  SUBROUTINE  RESETS  THE  FIVE  YIELD  VALUES  BASED 
UPON  THIS  ITERATION’S  MINIMUM  RESIDUAL. 
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**-K-******-»**')H(-***')H(")HHHH!-**»***-)l"tr  ■»■*■»**•»*******«**** 


COMMON  /ITEROT/  W <5> ,  P (5)  ,  DELTfiW (5) , DELTOP (5) , YLD <5) 

FIND  THE  MINIMUM  DELTA 
IF(DELTftW(5).LT.DELTflW<4) )  BO  TO  10 
IF(DELTflW(4).LT.  DELTAW<3) )  GO  TO  £0 
IF(DELTAN<3).LT. DELTAW{£) )  GO  TO  30 
IF(DELTAW(£) . LT.  DELTAW(l) >  GO  TO  40 

REDEFINE  YIELDS  BASED  UPON  THE  MINIMUM 

IF  DELTAW<1)  IS  MIN, 

DYLD=(W<£)-W(1> )*0.  £5 
W(S)=W(£) 

DELTAW (5) =DELTAW  <£) 

GO  TO  50 

IF  DELTAW (5)  IS  THE  MINIMUM, 

10  DYLD-(W<5)-W<4) )*0.  £5 
W(1)=W(4) 

DELTAW ( 1 ) =DELTAW ( 4 ) 

GO  TO  50 

IF  DELTAW (4)  IS  THE  MINIMUM, 

£0  DYLD=(W(5)-W<3.'i  >»0.  £5 
W(1)=W(3) 

DELTAW ( 1 ) =DELTAW  <3) 

GO  TO  50 

IF  DELTAW <3)  IS  THE  MINIMUM, 

30  DYLD=<W(4)-W(£) >*0. £5 
W(1)=W(£) 

W(5)=W(4) 

DELTAW  < 1 ) =DELTAW  <£> 

DELTAW (5) =DELTAW  <4) 

GO  TO  50 

IF  DELTAW (£)  IS  THE  MINIMUM, 

40  DYLD-<W<3)-W<1) )*0.  £5 
W<5)=W(3) 

DELTAW ( 5 ) =DELTAW ( 3 ) 

50  W<£)=W(l)+DYLb 
W(3)=W(£)+DYLD 
W(4)=W<3)+DYLD 
RETURN 
END 

SUBROUTINE  RESETP 
C  *****■»*•»*♦***■»******■»■)(■*♦*•«■********•»*******■»***» 
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C  THIS  SUBROUTINE  RESETS  THE  FIVE  PRESSURE  VALUES 

C  BASED  UPON  THIS  ITERATION’S  MINIMUM  RESIDUAL. 

C  **********************************-x-»**-innnnnnt** 

C 

C 

COMMON  /ITERAT/  W (5> , P (5) , DELTAW (5) ,  DELTAP <5) , YLD <5) 

r' 

u# 

C  FIND  THE  MINIMUM  DELTAP 

IF(DELTAP(5).LT. DELTAP<4) )  30  TO  10 
IF<DELTAP(4) .  LT.  DELTAPO)  )  SO  TO  £0 
IF (DELTAP (3). LT. DELTAP <2) )  GO  TO  30 
IF(DELTAP(£).LT.  DELTAP (1>>  GO  TO  40 
C 

C  REDEFINE  PRESSURES  BASED  UPON  THE  MINIMUM 

C 

C  IF  DELTAP (1)  IS  THE  MINIMUM, 

DPRESS=(P(£)-P(1) )*0.  £5 
P(5)=P(£> 

W(5)=W(£) 

DELTAP (5) =DELTAP (£) 

GO  TO  50 
C 

C  IF  DELTAP (S)  IS  THE  MINIMUM, 

10  DPRESS=<P<5)-P(4) )*.  £5 
P(1)=P(4) 

W<1)=W(4) 

DELTAP  < 1 ) =DELTAP (4) 

GO  TO  50 
C 

C  IF  DELTAP (4)  IS  THE  MINIMUM, 

£0  DPRESS=<P(5)-P<3) )  »0.  £5 
P(1)=P<3) 

W(1)=W(3) 

DELTAP ( 1 ) =DELT AP ( 3 ) 

GO  TO  50 
C 

C  IF  DELTAP (3)  IS  THE  MINIMUM, 

30  DPRESS= ( P ( 4 ) -P ( 2 ) ) *0.  £5 
P(1)=P(£) 

W(1)=W(2) 

DELTAP  ( 1  >  =-DELT AP  ( 2 ) 

P(5)=P(4) 

W(5)=W(4) 

DELTAP (5) =DELTAP (4) 

GO  TO  50 
C 

C  IF  DELTAP <2)  IS  THE  MINIMUM, 

40  DPRESS=(P(3)-P(1) )*0.  £5 
P<5)=P(3) 

W(5)=W(3) 
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c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 


c 

c 


c 

c 


DELTftP ( 5 ) =DELTflP ( 3 ) 

50  P<£)=P(1)+DPRESS 
P(3)=P<£>+DPRESS 
P(4)=P(3>+DPRESS 
RETURN 
END 

SUBROUTINE  RPNGE 

********************************************  t* 

THIS  SUBROUTINE  IS  AN  ITERATION  TO  FIND  THE  RANSE 
OF  THE  ESTIMATED  PEAK  PRESSURE  FDR  THE  ESTIMATED 
YIELD.  THIS  IS  NECESSARY  FOR  COMPUTATION  OF  THE 
SPEICHER-BRODE  PRESSURE  HISTORY,  TIME  OF  ARRIVAL 
AND  POSITIVE  PHASE  DURATION. 

**»**«-»****************«*************it*********** 


COMMON  /ESTIM  /  PSOI, WI, PP, W13, PSOF, WF,  FSB 
COMMON  /PEAK  /  DP,  TA, PSO, ALPF 
COMMON  /SBCONS/  RSKFT, YS, S, XM 
COMMON  /COUNT  /  ICOUNT, lOPT, LFILT 

INITIAL  RANGE  SPREAD 
IFdOPT.  NE.  3)  GO  TO  78 
PP=PSOI 

W13=WI*». 33333 
78  R1=0. 01 
R£=0. 1 
R3=l. 0 
RA=10. 

HOB  EQUAL  TO  ZERO 
Y=0. 

YS1=0. 

YS£=0. 

YS3=0. 

YS4=0. 

DO  100  1  =  1,  1000 
RS1=R1/W13 
RS£=R£/W13 
RS3=R3/W13 
RS4=R4/W13 

CALCULATE  PSO  FOR  EACH  TRIAL  SCALED  RANGE 
CALL  PPEAK(RS1, YSl, PI) 

DP1=DP 

CALL  PPEAK<RS£, YS£,  P£) 

DP£=DP 

CALL  PPEAK ( RS3, YS3, P3) 

DP3=DP 

CALL  PPEAK ( RS4, YS4, P4> 


on  no  on  on 


DPA=DP 


FIND  BOUNDING  RfiNGES 

IF(PP.  GT. P2  .AND.  PP.  LT.  PI)  GO  TO  110 
IF<PP.  6T.P3  .AND.  PP.  LT.  P£)  GO  TO  120 
IF(PP. GT.  P4  .AND.  PP.  LT. P3)  GO  TO  130 
WRITE(6, 11A0) 

11A0  FORMAT (2X, *PRESSURE  OUT  OF  RANGE*) 

STOP 11 

BETWEEN  R1  AND  R£ 

110  DR=<R2-Rl)/3. 

RA=R2 
R£=R1+DR 
R3=R2+DR 
GO  TO  99 

BETWEEN  R£  AND  R3 
120  DR=(R3-R£)/3. 

R1=R2 
RA=R3 
R2=R1+DR 
R3=R2+DR 
GO  TO  99 

BETWEEN  R3  AND  RA 
130  DR=(RA-R3)/3. 

R1=R3 

R£=R14-DR 

R3=R2+DR 

99  IF( (RA-R1).LE..0001)  GO  TO  101 

100  CONTINUE 
WRITE (6, 1100) 

1100  FORMAT (2X,  *FAILED  TO  CONVERGE  ON  RANGE*) 

WRITE <6, 1200)  I 

1200  FORMAT <£X, *1=*, 15) 

WRITE(6,  1201)  PP,  R1,RA 

1201  FORMAT  <2X, *PP=*,  E12.  5,  /,  £X,  *R1=*,  E12.  5,  /,  2X,  *RA=*,  E12. 5) 
STOP 12 

101  RAKFT=(R1+R£:+R3+RA)*0.  £5 
RSKFT=RAKFT/W13 

DP=  <DP1+DP2+DP3+DPA) *0. £5 
FS0=1. /<DP/1000. ) 

IFdCOUNT.NE.  1)  GO  TO  103 
TASEC=(TA/1000.  )*W13 
DPOS=(DP/1000.  )*W13 
RANKM=RAKFT*.  30Aa 
PSOF=PP 

WP=W13*W13*W13 

C 
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WRITE  FINfiL  RESULTS  TO  OUTPUT  FILE 
WRITE (6, 1 10£) PSOF, WF, RPNKM,  TftSEC,  DPOS 
,10£  FORMftT(//, IX, *++++++++++++++++++++++++++++++++++h+++++++*, /, 

*  £X,*PEftK  OVERPRESSURE,  MPft==*,6X,El£.  5, //, 

*  £X,*NUCLEfiR  YIELD,  KT=*,  1 IX,  E1£.S, //, 

*  £X,  GRANGE  FROM  GZ,  KM=*,  1 1 X,  El£.  5, //, 

*  £X,*TIME  OF  ARRIVAL,  SEC==*,8X,E1£.  5, //, 

*  £X,»POSITIVE  PHASE  DURATION, SEC=*, El£. 5, /,  1 X, 

WRITE (48, 1 103)  PSOF, WF 
WRITE (48, 1104)  DP, TA, RSKFT 
.103  FORMAT  (£E15.  8) 

.104  FORMAT (3E15. 8) 

103  RETURN 
END 

SUBROUTINE  PPEAK (X, Y, PEAKP) 

THIS  SUBROUTINE  CALCULATES  THE  PEAK  OVERPRESSURE  (MPA), 

TIME  OF  ARRIVAL  (TA, MS/KT**l/3) ,  AND  POSITIVE  PHASE 
DURATION  (DP, MS/KT**l/3)  AFTER  SPEICHER-BRODE,  JUNE,  198£. 

****•»*«■«••»**************•»  X  **•»*****■»*■♦*•)(•*•»■»****•»•******•»*■»*** 


COMMON  /PEAK  /  DP, TA, PSO, ALPF 
COMMON  /SBCONS/  RSKFT, YS, S, XM 

XLEAST=1. E-g 
YLEAST=1. E-9 
ZMAX=100. 

IF(X. LT. XLEAST)  X=XLEAST 
IF(Y. LT. YLEAST)  Y=YLEAST 
R=SQRT(X*X+Y*Y) 

R£=R*R 

R3=R*R£ 

R4=R£*R£ 

R6=R£*R4 

R8=R4*R4 

Z=Y/X 

Z£=Z*Z 

Z3=Z*Z£ 

Z5*Z£*Z3 

Z17=*Z**17. 

Z18=Z**18. 

Y7=Y**7. 

IF(Z.GT. ZMAX)  Z=ZMAX 

XM=170. *Y/ (1. +337. *Y**. £5) +.  914*Y**£.  5 


SCALED  TIME  OF  ARRIVAL 


n  n  n  non 


Ul=(. 543-21. 8*R+386.  *R2+E383.  *R3)*R8 
U2=2. 99E-14-1. 91E-10*R2+1. 03EE-6*R4-4.  43E-6*R6 
U3= ( 1 . 028+2. 087*R+2.  69*R2) *R8 
UTft=Ul/ (U2+U3) 

TA=UTft 

IF<X. LT. XM)  GO  TO  101 

Wl=  < 1 . 086—34. 605*R+486. 3*R2+E383. *R3) *R8 
W2=3.  0137E-13-1.2128E-9*R2+4.  128E-6-»frR4-l .  11&E-S*R6 
W3=<1. 632+2. 629*R+2.  69*R2)*R8 
WTft=Wl/ (W2+W3) 

TA=LJTAKXM/X+WTft*(l.  -XM/X) 

SCALED  POSITIVE  PHASE  DURATION 

101  S=l.-1. 1E10*Y7/(1.+1.  lE10*Y7)-<2.  441E-8*Y*Y/ 

*  (1. +9. E10*Y7) )*(1.  / (4.  41E-11+X**10.  ) ) 

DP= ( ( 1640700. +24629. *TA+41&.  15*TA*TA) / 

*  ( 1 0880. +619. 76*T A+T A*T A ) ) 

*  *<. 4+. 001204* (TA**1. 5) /<!.+. 001559*TA**1 . 5)  + 

*  (.0426+. 5486*  <TA**.  25) / (1.  +.  00357*TA**1 . 5) ) *S) 

AA=1. 22- (3. 908*Z2) / < 1.  +810.  2*Z5) 

BB=2. 321  + (Z 18/ <1.  +1.  113*Z18) )*6.  195- ( . 03831*Z 17) / 

*  (l.+.  02415*Z17)+.  6692/(1. +4164.  *Z**a.  > 

CC^4. 153- (1.  149*Z18)/(1.+1.641*Z18)-1.  1 / ( 1 . +2.  771*Z**2. 5) 
DD=-4. 166+ (25. 76*Z**1. 75) / ( 1 .  +1 . 382*Z 18) +8.  257*Z/ (1. +3. 219*Z) 
EE= 1 . - ( . 004642*Z 1 8 ) / ( 1 .  +.  003886*Z 1 8 ) 

FF=. 6096+ (2.  879*Z**9.  25) / ( 1.  +2.  359*Z**14.  5) -17.  15*Z2/ 

*  (1. +71. 66*Z3) 

GG=1. 83+5. 361*Z2/ ( 1.  +.  3139*Z**6.  ) 

HH=-(64. 67*Z5+. 2905) / (1.  +441 . 5*Z5) -1 .  389*Z/ (1. +49. 03*Z5)  + 

*  (8. 808*Z**1. 5) / (1.  +154.  5*Z**3.  5)  +  ( .  0014*R2/ ( 1 .  -. 158*R+ 

*  . 0486*R**1. 5+.  00128*R2) )*(1.  /(1.+2. *Y)) 

PEAK  OVERPRESSURE 

P0=10.  47/ ( R**AA) +BB/ ( R**CC) +DD*EE/ ( 1 .  +FF*R**GG) +HH 
PEAKP=P0*.  006894757 
RETURN 
END 

SUBROUTINE  SPBRODE 
C  *■)(■***#**#**•»***********#****************************** 

C  THIS  SUBROUTINE  CALCULATES  THE  PRESSURE  HISTORY  FOR 

C  THE  FINAL  PRESSURE-YIELD  PAIR  DETERMINED  BY  SUBROUTINE 

C  FIT.  IT  USES  THE  SPEICHER-BRODE  JUNE,  1982  ALGORITHM. 

C  *****■)(•#****■«■*****■)»•**•»•»*********•»********************** 

C 

C 

COMMON  /THIST  /  TTIM (6000) ,  PRESS ( 12000)  ,  TIMP (2999) , PIMP (2999) , 

*  PFILT(6000) 

COMMON  /FFT  /  FRQ ( 3001 ),  AMP (3001 ),  XFFT (3001 ) 
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c 


c 

c 

c 


COMMON  ZEST I M  / 
COMMON  /PEAK  / 
COMMON  ZFILT  / 
COMMON  ZSBCONSZ 
COMMON  /POINTS/ 
COMMON  /IMP  / 
COMMON  /COUNT  / 
COMMON  /PLOTV  / 
COMPLEX  XFFT 
DIMENSION  IWKB( 
DATA  JCOUNT/0/ 


PSOI,  WI,  PP, W13, PSOF, WF, FS0 
DP,  TA,  PSO,  ALPF 

IFILT, FL0(7) , PFDMX  <7) , PFBMX  <7) 

RSKFT,  YS,  S,  XM 

NEPTS,  NBPTS,  NI, NEF, NBF 

I  IMP,  DTD,  DTB, TPEB, DTBN 

I COUNT,  I OPT, LFILT 

ITL(8) , ISTL(a) , IDB 

11) 


IF<IOPT.NE.  3)  SO  TO  5 
I TL ( 1 ) = 1 0HC ALCUL ATED 
ITL(£)=10H  SPEICHER 
ITL(3)=10HBRODE  PRES 
ITL(4)=10HSURE  HISTO 
ITL(5)=10HRY 
ITL (6) =10H 
ITL(7)=10H 
ITL(a)-10H 

WRITE<4a,£6)  <ITL(IO),  IO=l,a) 

£6  FORMAT (aA10) 

CALCULATE  SPEICHER-BRODE  TIMESTEP  BASED 
UPON  THE  POSITIVE  PHASE  DURATION. 
DTB=DP/NBPTS 


GO  TO  15 

5  ISTLd  )=10HWITH  FOURF 
ISTL(£)=10HIT  SPEICHE 
ISTL(3)=10HR  ERODE 
ISTL<4)=10H 
ISTL(5)=10H 
ISTL(&)=10H 
ISTL<7)=10H 
ISTL<8)=10H 

WRITE(A8, £&)  (ISTL(IG), 1G=1,8) 

CALL  FMAX (PRESS, NEPTS, YPMN, YPMX) 

CALL  FMAX (TTIM, NEPTS, XPMN, XPMX) 

WRITE (48, £00)  NEPTS, XPMN, XPMX, YPMN, YPMX 
WRITE(4a, £10)  (TTIM(IU), IU=1,  NEPTS) 
WRITE (48, £10)  (PRESS (IP), IP=1, NEPTS) 

£00  F0RMAT(I5,4E15.a) 

£10  FORMAT (10E15. 8) 

ICOUNT=0 

C 

C  FIND  THE  PEAKS  OF  THE  LOW  PASS 

C  FILTERED  DATA  PRESSURE  HISTORIES 

DO  7  1  =  1,7 

CALL  FILTER (DTD, NEPTS) 


35 


non  nnnnnn 


CALL  FMAX (PFILT,  NEPTS,  PFDMN,  PFDMX ( I ) ) 

7  CONTINUE 
ICOUNT=l 

CALCULATE  SPEICHER-BRODE  TIME  STEP  BASED 
UPON  THE  DATA  TIME  STEP  FOR  FILTERINB 
DTB=DTD*1000. /W13 

CALCULATE  THE  SPEICHER-BRODE  TIME  STEP  BASED 
UPON  THE  POSITIVE  PHASE  DURATION  FOR  OVERLAYS 
35  IF(JCOUNTc EQ. 1)  DTB=DP/NBPTS 
15  DO  £5  KJ=1,NBPTS 
TTIM(KJ)=0. 

PRESS (KJ)=0. 

£5  CONTINUE 
X=RSKFT 
TF=TA+DP 
P0=PSOF*145. 038 

F- < . 01477*  <TA**.  75) /(!.+.  00583S*TA) +7.  40£E-5* (TA**£. 5) / 

*  < 1 . +1. 4£9E-a*TA**4.  75) £16) *S+.  7076-3.  077E-5* 

*  TA*TA*TA/ ( 1 .  +4.  367E-5*TA*TA*TA) 

G= 1 0.  + ( 77. 58-64.  99* (TA**.  1 £5 ) / ( 1 .  +.  0434a*SQRT  < TA ) ) ) *S 
H=£. 753+. 05601 *TA/ <1.  +1 . 473E-9*TA**5.  )  +  (.  01769*TA/ 

*  < 1. +3. £07E-10*TA**4.  £5) -.  03£09* (TA**1.  £5) / < 1.  +9.  914E-a* 

*  TA**4. )-l. 6) *S 

CALCULATE  PRESSURE  HISTORY 

DO  400  J=1,NBPTS 

T=TA+(J-£)*DTB 
C  SAVE  UNSCALED  TIMES 

TT I M ( J ) =T*W 13/1 000. 

PRESS <J)=0. 

IF<T.LT.  TA)  SC  TO  400 
IF<T.GT.TF)  GO  TO  410 

B=  <F*  <TA/T) **G+ ( 1.  -F) * (TA/T) **H) * ( 1.  - <T-TA) /DP) 

XE=3. 039*Y/ (1.  +6. 7*Y) 

E=ABS< (X-XM) / <XE-XM) ) 

IF(E.  GT.  50.  )  E=50. 

D=. £3+583000.  *Y*Y/ <£6667.  +1.  E6*Y*Y)+.  £7*E+ (.  5-583000.  *Y*Y/ 

*  (£6667. +1. E6*Y*Y) )*E**5. 

A=  <D-1. ) * (1.  - <E**£0.  ) / ( 1.  +  (E**£0.  ) ) ) 

DT=474.£*Y*(X-XM)**1. £5 
IF(DT.  LT.  1.  E-9)  DT=l.E-9 
GA=(T-TA) /DT 
IFCGA. GT. 400. )  GA=400. 

V=l.  +  (3. £8E11* (Y**6.  ) / (1.  +1. 5E1£*Y**G.  75) )*(GA*GA*GA/ 

*  (6.  13+GA*GA*GA) )* ( 1.  / ( 1.  +9.  £3*E*E) ) 

C=( (1. 04-£40.  9*(X**4) / (1.  +£31. 7*X**4) )* (GA**7) / 

*  <(!.+.  9£3*GA**8.  5)*<l.+A)))*(l.-(  (T-TA)  /DP)**8.  ) 
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*  *£. 3E13*Y**9.  / <1- +£.  3E13*Y**9) 

POFT=P0* ( 1 .  +R) * (B*V+C) 

IF(X.LT.  XM.  OR  .Y.GT.  .38)  POFT=P0*B 
PRESS (J)=PaFT/ 145. 

400  CONTINUE 
C 

410  JC0UNT=JC0UNT+1 

UNSCRLE  THE  SPEICHER-BRODE  TIMESTEP 
DTBN=DTB*W 13/1 000. 

IF(JCOUNT.GT.  1  .OR.  I  OPT.  ED.  3)  GO  TO  900 

FIND  THE  PERKS  OF  THE  LOW  PASS  FILTERED 
SPEICHER-BRODE  PRESSURE  HISTORIES 
LFILT=0 
DO  17  J=l,7 

CALL  FILTER (DTBN,  NBPTS) 

CALL  FMAX (PFILT, NBPTS, PFBMN, PFBMX ( J) > 

17  CONTINUE 

FIND  THE  LOW  PASS  FIDELITY  FREQUENCY 

DO  £7  K=l,7 

PFMAX=PFDMX  <K) *0.  90 
IF <PFMAX. LE. PFBMX (K) )  GO  TO  47 
£7  CONTINUE 
WRITE <6, 37) 

37  FORMAT (£X, *+++  FAILED  TO  LOCATE  LOW  PASS  FIDELITY  +++*) 
ALPF=-999. 

WRITE (48, 57)  ALPF 
GO  TO  35 
47  ALPF=FLO<K) 

WRITE <48, 57)  ALPF 
57  FORMAT (F10. 0) 

WRITE (6, 67)  ALPF 

67  FORMAT <£X, *+++  LOW  PASS  FIDELITY  (HZ)  =  *, F10. 0, *  +++*) 
IF<JCOUNT.Ea. 1)  GO  TO  35 

DETERMINE  NUMBER  OF  SPEICHER  BRODE  PAIRS  TO 
BE  PLOTTED  FOR  OVERLAY 
900  TE=NEPTS*DTD 

NPPTS=IFIX (TE/DTBN) 

IF ( I OPT. EQ. 3)  NPPTS=NBPTS 
WRITE (48, 450)  NPPTS 
450  FORMAT (15) 

IFdOPT.  EQ.  3)  GO  TO  310 

AFFECT  A  TIME  SHIFT  IN  SPEICHER-BRODE  HISTORY 
TO  ALLOW  THE  OVERLAY  TO  BE  PROPERLY  PHASED 
TSHFT=(TA*W13/1000. ) -TPEB 


37 


non  nnnll  nn 


DO  800  JT=1,NBPTS 

TT I M ( JT ) =TT I M ( JT ) -TSHFT 
800  CONTINUE 

BO  TO  130 

810  CftLL  FMftX (TTIM,  NBPTS,  XPMN, XPMX) 

COLL  FMflX (PRESS,  NBPTS,  YPMN, YPMX) 

WRITE (48, 840)  XPMN,  XPMX,  YPMN, YPMX 
840  FORMAT (4E 15. 8) 

IFdFILT.  LT.  0)  GO  TO  130 

CALL  FOR  FILTERS  TO  BE  EXECUTED 
CALL  FLOOP ( TT I M,  PRESS,  DTBN, NBPTS, PF I LT ) 

RETURN 

130  WRITE(48, £10)  (TTIM ( IJ) ,  IJ=1 , NPPTS) 

_  WRITE(48.£10)  (PRESS(JI) .  JI  =  1. NPPTS) 

WRI FE (49, £00)  NPPTS,  DTBN 
WRITE(49, £10)  (TTIM(IJ),  IJ=1, NPPTS) 

WRITE (49, £10)  (PRESS (JI ) ,  JI  =  1. NPPTS) 

IFHOP'T.  E(J.  3)  GO  TO  850 

IMPULSE 

135  ITL(5)=10H  IMPULSE  H 
ITL(6)=10HISTORY 
WRITE(4a, £15)  ITL(5) , ITL(6) 

£15  FORMAT (£A10) 

CALL  FMAX(TIMP, NI,XIMN, XIMX) 

CALL  FMAX (PIMP, NI, YIMN, YIMX) 

WRITE (48, £00)  NI,  X IMN,  X IMX, YIMN, YIMX 
WRITE(48,£10)  (TIMP(IY), IY=1,NI) 

WRITE(4a,£i0)  (PIMP(IT),  IT=1,NI) 

850  IIMP=£ 

CALL  IMPULSE ( I  IMP,  DTBN,  NPPTS, NI ) 

WRITE (48, 450)  NI 
IF  (lOPT. NE. 3)  GO  TO  150 
ITL(3)=10HBRODE  IMPU 
ITL(4)=10HLSE  HISTOR 
ITL(5)=10HY 

WRITE (48, ££5)  ITL(3), ITL(4), ITL(5) 

£25  FORMAT (3A10) 

CALL  FMAX (TIMP, NI, XIMN, XIMX) 

CALL  FMAX (PIMP, NI, YIMN, YIMX) 

WRITE(48, 840)  XIMN,  XIMX,  YIMN, YIMX 
150  WRITE(48, £10)  (TIMP (KJ) ,  KJ=1 , NI > 

WRITE(48,£10)  (PIMP(KL) , KL=1, NI) 

IFdOPT.  NE.  1)  GO  TO  175 

FIND  THE  FOURIER  TRANSFORM  AND  CALCULATE  AMPLITUDE. 
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ITL<5)=10H  FOURIER  ft 
ITL(6)=10HMPLITUDE  S 
ITL(7)=10HPECTRUM 

WRITE(Aa, 2£5)  ITL<5), ITL(6>, ITL(7) 

CftLL  FMfiX <FRQ, NEF, XFMN,  XFMX) 

CftLL  FMflX (AMP,  NEF, YFMN,  YFMX) 

WRITE  (48,  200)  NEF,  XFMN,  XFMX,  YFMN, FMX 
WRITE<4a,  210)  (FRQdO),  10=1,  NEF) 

WRITE(4a,  210)  (AMP(IP) ,  IP=1, NEF) 

175  TOTT=DTBN*NBPTS 
C  FREQUENCY  INCREMENT 

DFB=1. /TOTT 
FQB=0. 

WKB=0. 

NBF=NBPTS/2+l 
DO  349  LK=1,NBF 
FRQ(LK)=0. 

AMP(LK)=0. 

XFFT(LK)=0. 

349  CONTINUE 

CALL  FFTRC (PRESS,  NBPTS,  XFFT,  IWKB, WKB) 

C  AMPLITUDE  SPECTRUM 

DO  500  KK=1,NBF 
FQB=FQB+DFB 
FRQ(KK)=FQB 

XRB=REAL(XFFT(KK) ) /NBPTS 
XIB=AIMAG(XFFT(KK) ) /NBPTS 
AMP (KK) =SQRT (XRB*XRB+XIB*X IB) *TOTT 
500  CONTINUE 
C 

WRITE (48,  450)  NBF 
IFdOPT.  NE.  3)  BO  TO  165 
ITL(3)=10HBRODE  FOUR 
ITL(4)=10HIER  AMPLIT 
ITL(5)=10HUDE  SPECTR 
ITL(6)=10HUM 

WRITE(4a,  £35)  ITL(3),  ITL(4), ITL(5),  ITL(6) 

£35  FORMAT (4A10) 

CALL  FMAX (FRQ,  NBF, XFMN, XFMX) 

CALL  FMAX (AMP,  NBF,  YFMN, YFMX) 

WRITE (48, 840)  XFMN,  XFMX,  YFMN, YFMX 
165  WRITE(48,  £10)  (FRQ  dU) ,  IU=1,  NBF) 

WRITE(48,  210)  (AMPdE),  IE=1,  NBF) 

RETURN 
END 

SUBROUTINE  FLOOP (TTIM,  PRESS, DT, NP, PFILT) 
****•»*******■»***************•************•»***•)(■■» 
THIS  SUBROUTINE  PERFORMS  THE  LOOPING  REQUIRED 
TO  FILTER  THE  DATA  OR  THE  BRODE  UP  TO  SEVEN 
TIMES.  FOR  LESS  THAN  SEVEN  FILTER  LEVELS, 


39 


nn  n  nnnnn 


FLO  MUST  BE  SET  TO  0.  IN  THE  INPUT  DECK  IN 
ORDER  TO  ESCAPE  THE  LOOP. 

*«********«***************-inf-****************** 


COMMON  /FILT/  IFILT,  FLO (7) ,  PFDMX (7) , PFEMX <7) 
DIMENSION  TTIM  < 1 ) , PRESS ( 1 ) , PFILT ( 1 ) 

DO  750  JF=1,7 

IF(FLO<JF).EG!.  0.  )  GO  TO  555 
IFLAG=1 

WRITE <48, 95)  IFLAG 

95  FORMAT (15) 

WRITE <48, 96)  FLO<JF) 

96  FORMAT <F10.0) 

DO  7£5  KF=1,NP 

PFILT (KF)=0. 

7S5  CONTINUE 


CALL  TO  FILTER 

CALL  FILTER <DT,  NP) 

CALL  FMAX  <PFILT,  NP,  YFMN, YFMX ) 
WRITE <48, 100)  YFMN,  YFMX 
100  FORMAT <£E1S. 8) 

WRITE<48, 105)  <TTIM<LF) , LF=1, NP) 
WRITE  <48, 105)  <PFILT  <MF) , MF=1, NP) 
105  FORMAT <10E13.  8) 

750  CONTINUE 
555  IFLAG=-i 

WRITE <48, 95)  IFLAG 

RETURN 

END 


SUBROUTINE  SPLINE  <TLAST,  NP,  TTIM, PRESS) 

C  ****»************«****#**'»***-«"»-*4f***********»** 

C  THIS  SUBROUTINE  SETS  UP  A  COSINE  SQUARED  SPLINE 

C  FUNCTION  AND  APPLIES  IT  TO  THE  FINAL  15X  OF  THE 

C  PRESSURE  HISTORY  TO  AVOID  A  FREQUENCY  IMPULSE 

C  IN  TRUNCATED  RECORDS. 

C  **********X"»*-»****-)HHH<"»*************)HHH{-*****-»* 

C 

DIMENSION  TTIM<1),  PRESS<1) 

C 

PIE=3. 14159£7 
K=IFIX<.85*NP) 

N=NP-K+1 

T1=TTIM<K) 

DO  10  J=1,N 

TFACT=  <TTIM  <K) -T1 ) / <TLAST-T1 ) 

SFACT=COS  <TFACT*PIE*. 5) 

SFACT=SFACT*SFACT 
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c 

c 

n 

w 

c 

c 

c 

c 


c 


c 

c 

c 

c 

c 

c 

c 

c 


PRESS  <  K ) =PRESS  <  K  >  *SFACT 
K=:K+1 

10  CONTINUE 
RETURN 
END 

SUBROUTINE  IMPULSE ( I IMP, DT, NP, NI) 
******•****•»•»*#*******♦*******■!«■****#*■»***#♦**•!«■******•»* 
THIS  SUBROUTINE  CftLCULPTES  THE  IMPULSE  OF  THE  INPUT 
PRESSURE  DATfi  (IIMP=1)  OR  OF  THE  CPLCULftTED  SPEICHER- 
BRODE  <IIMP=£)  BY  SIMPSON’S  APPROX  I MOTION. 


COMMON  /THIST  /  TTIM <6000) ,  PRESS ( 12000) ,  TIMP (£999) ,  PIMP <£99S) , 
*  PFILT<6000) 


NTMP=NP-3 

NI=NTMP/£ 

DO  90  1=1, NI 

TIMP(I)=0. 

PIMP(I)=0. 

90  CONTINUE 
IJ=0 

SUMIMP=0. 

DO  80  J=3,  NTMP,  £ 

IJ=IJ+1 

T1MP<IJ>=TTIM<J) 

AREA=  <PRESS ( J-1 ) +4. *PRESS ( J) +PRESS  < J+1 ) ) *DT/3. 
SUM I MP=SUM I MP+ A  RE A 
PIMP<IJ)=SUMIMP 
80  CONTINUE 
RETURN 
END 

SUBROUTINE  FILTER <DT, NP) 

THIS  SUBROUTINE  FILTERS  THE  INPUT  PRESSURE  HISTORY 
(DATA  OR  SPEICHER-BRODE).  IT  USES  THE  DIFFERENCE 
EQUATIONS  DERIVED  FOR  A  SECOND  ORDER  BUTTERWORTH 
FILTER  AS  PRESENTED  BY  STEARNS,  1975. 


COMMON  /THIST 

COMMON  /COUNT 
COMMON  /FILT 
DATA  LFILT/0/ 
Pl=3. 1415927 
S£=SQRT(£.  ) 
LFILT=LFIlT+1 


/  TTIM (6000) , PRESS ( 12000) , TIMP (£999) 
■pFILT<6000) 

/  ICOUNT,  lOPT, LFILT 
/  IFILT, FLO (7) , PFDMX (7) , PFBMX (7) 


PIMP(£999), 
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LOW  PftSS  FILTER  COEFFICIENTS 

ftT=T AN ( P I *FLO  <  LF I LT ) *DT ) 

ATe=AT*AT 
Al  =  la  +S£*AT+AT£ 

A=AT£/A1 
Bl=£.  *(AT£-1.  > 

B=B1/A1 

Cl=l. -S£*AT+AT£ 

C=C1/A1 
FAC=1. 

CALCULATE  THE  FILTERED  HISTORY 
150  PFILT(1)=A*PRESS<1) 

PFILT (£) =A* (PRESS (£) +£*FAC*PRESS ( 1 ) ) -B*PFILT ( 1 ) 

DO  £00  1=3, NP 

pC=A^<PRESS(I)+£.  *FAC*PRESS(I-l)+PRESS(I-£) ) 
PF I LT ( I ) =PC-B*PF I LT ( I - 1 ) -C*PF I LT ( I -£ ) 

£00  CONTINUE 
RETURN 
END 

SUBROUTINE  FMAX ( ARY, NA,  XMN, XMX ) 

C  »***■»*#*■)(■*♦*•«•♦*■****■»**■»*********■«■**■«■**#**•»■)(•***** 

C  THIS  SUBROUTINE  FINDS  THE  MAXIMUMS  AND  MINIMUMS 

C  OF  THE  VARIOUS  ARRAYS  TO  BE  PLOTTED  BY  FOURPLT 

C  *■»*****■***♦**■»*•»■»***•»*■»*•»**«****»***■»******■*•***■» 

C 

c 

DIMENSION  ARY(NA) 

C 

XMN  =  ARY(l) 

XMX  =  ARY(l) 

IF(NA. EQ. 1)  RETURN 
DO  10  I=£,NA 

IF(XMN.GT. ARY(I) )  XMN  =  ARY ( I ) 

IF<XMX.LT. ARY(I) )  XMX  =  ARY ( I ) 

10  CONTINUE 
C 

RETURN 

END 
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APPENDIX  B 


PROGRAM  FREQRES  USER'S  MANUAL 


B.l  INPUT  VARIABLES 

Table  2  lists  all  variables  used  in  the  input  file  for  program 
FREQRES  and  the  format  in  which  they  occur.  In  order  to  run  FREQRES, 
first  one  has  to  select  two  digitized  time  history  records  with  the  same 
time  step  and  record  duration.  One  of  the  records  should  be  a  loading 
waveform  and  is  considered  the  input.  The  other  record  should  be  some 
type  of  response  time  history  and  is  considered  the  output.  From 
Table  2,  NEPTS  is  the  number  of  discrete  data  points  in  each  of  the 
digitized  records.  NSKIP  is  the  skip  factor  if  the  user  wishes  to  work 
with  less  than  NEPTS  points.  The  total  number  of  points  from  the  input 
and  output  data  records  now  considered  for  analysis  is 
NPT=NEPTS/NSKIP 

The  digitized  time  step  is  multiplied  by  NSKIP  to  get  a  new  resultant  time 
step.  The  discrete  data  points  are  averaged  locally,  so  that  all  NEPTS 
points  read  from  tape  are  considered  in  the  analysis.  TFAC,  XFAC,  and 
AFAC  are  the  time,  input  data,  and  output  data  conversion  factors,  if  the 
user  does  not  want  to  work  .■  units  as  specified  on  the  data  record 
tapes.  Also,  a  -1.0  for  XFAC  or  AFAC  will  invert  the  input  or  output  data 
records,  respectively. 

ISPBX  is  a  trigger  indicating  whether  or  not  the  user  wants  to  spline 
the  beginning  portion  of  the  input  time  history  with  a  cosine  squared 
spline,  such  that  the  input  value  at  zero  time  is  forced  to  zero.  If  one 
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elects  to  apply  this  spline,  then  TSPBX  is  the  time  before  which  the 
spline  is  applied  to  the  input  data.  ISPEX  is  a  trigger  indicating 
whether  or  not  the  user  wants  to  apply  the  cosine  squared  spline  to  the 
final  portion  of  the  input  time  history  such  that  the  final  input  value  is 
forced  to  zero.  TSPEX  specifies  the  time  at  which  the  final  portion 
spline  begins.  If  a  final  portion  spline  is  requested  and  TSPEX  is  left 
blank,  then  the  default  is  the  final  15  percent  of  the  input  record 
duration  is  splined.  ISPBA,  TSPBA,  ISPEA,  and  TSPEA  describe  the  spline 
conditions  for  the  output  data  record. 

From  card  3  of  Table  2,  IBLX  is  the  baseline  correction  trigger  for 
the  input  data.  Either  a  constant  or  a  linear  baseline  correction  can  be 
applied.  DELPX  specifies  the  amount  of  baseline  correction.  SBX  is  the 
start  time  for  the  baseline  correction  and  EBX  is  the  end  time.  If  SBX 
and  EBX  are  specified  as  equal,  then  the  full  value  of  DELPX  is  added  to 
all  input  values  after  time  SBX.  If  EBX  is  greater  than  SBX,  then  the 
input  plot  is  rotated  about  the  point  defined  at  SBX  by  the  amount  DELPX 
at  time  EBX.  The  resulting  input  data  trace  can  then  be  used  as  is  or 
integrated  according  to  the  integration  trigger  INTX.  IBLA,  DELPA,  SBA, 
EBA,  and  INTA  describe  the  baseline  correction  and  integration  conditions 
for  the  output  data. 

Cards  5  through  16  in  Table  2  specify  labels  for  the  X  and  Y  axis 
for  all  of  the  different  types  of  plots  resulting  from  a  FREQRES  run. 

Table  2  describes  the  uses  of  each  label.  All  of  the  labels  should  be 
centered  within  the  first  30  columns  of  each  line  of  the  input  file  so 
that  the  labels  are  centered  on  the  axis  of  the  plots. 
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B.2  PROGRAM  STRUCTURE 

A  listing  of  program  FREQRES  is  provided  in  Appendix  C.  Figure  33 
provides  a  flow  chart  of  the  program  structure.  It  is  easy  to  follow  the 
program  listing  as  one  preceeds  through  the  flow  chart,  since  the  program 
Is  well  documented  with  comment  cards.  The  first  part  of  the  program  has 
extensive  comment  cards  describing  the  program  input  file.  The  first 
executable  part  of  the  program  reads  the  entire  input  file  and  then 
documents  all  input  file  conditions  in  the  printed  output  file.  The 
program  then  writes  the  input  file  conditions  to  the  time  history  plot 
file.  The  input  data  record  time  history  is  read  from  TAPE26  and  then  can 
be  baseline  corrected,  integrated,  and/or  splined  according  to  input  file 
specifications.  The  resulting  input  data  record  time  history  is  then 
stored  on  a  plot  file  (TAPE48)  for  post-process  plotting.  The  output  data^ 
record  time  history  is  read  from  TAPE27  and  it  can  also  be  baseline 
corrected,  integrated,  and/or  splined  according  to  input  file 
specification.  The  resulting  output  data  record  time  history  is  then 
stored  on  the  plot  file.  Another  input  data  record  of  arbi tarary 
specification  (such  as  a  "best-fit"  Speicher-Brode  waveform  to  a  HEST 
pressure  input  data  record)  is  read  from  TAPES.  If  the  time  step  from  the 
new  input  data  record  is  not  within  a  1  percent  tolerance  level  of  the 
time  step  for  the  input  and  output  data  records  read  earlier,  then  the 
program  will  stop.  The  time  steps  should  be  nearly  identical  for  best 
results  in  this  analysis  procedure.  As  an  example,  assume  that  the  time 
step  for  input  and  output  data  records  is  5  x  10"®  seconds  and  the  time 

A 

Step  for  the  new  input  data  record  is  1.337  x  10"“^  seconds.  Then  one 
would  have  to  specify  an  NSKIP  value  of  (1.337  x  10"^)/(5  x  10"®^  =  27 
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in  order  to  allow  the  program  to  run.  The  resulting  time  steps  would  be 
1.35  X  10”^  seconds  and  1.337  x  10"^  seconds,  which  are  within  the  1 
percent  tolerance  level.  If  the  1  percent  tolerance  level  is  too 
restrictive  for  some  data  record  combinations,  the  user  may  have  to  relax 
the  tolerance  level  to  2  percent,  or  at  the  most  3  percent,  by  updating 
the  program. 

The  program  FFT's  the  input  and  output  data  record  time  histories 
(see  Section  2.2)  and  stores  the  FFT  amplitude  spectrums  in  the  plot 
file.  The  FRF  is  calculated  by  dividing  the  output  record  FFT  by  the 
input  record  FFT  (see  Section  2.4).  The  FRF  phase  angle  information  is  of 
little  interest  in  this  analysis,  but  FRF  amplitude  ratios  are. 

Therefore,  the  FRF  amplitude  ratios  are  saved  in  the  plot  file. 

As  was  discussed  in  Section  2.3,  the  inverse  FFT  of  an  FFT  does  not 
give  back  the  same  exact  discrete  time  series.  An  inverse  fft  is  applied 
directly  to  the  output  record  FFT,  in  order  to  get  an  output  record  time 
history  which  includes  these  alterations. 

The  program  then  FFT's  the  new  input  data  record  read  earlier  from 
TAPE3.  This  FFT  is  then  multiplied  by  the  FRF  through  complex  math  to 
obtain  a  new  modified  output  response  FFT.  The  program  then  inverse  FFT's 
this  new  FFT  to  obtain  a  modified  output  response  time  history. 

B-3  SAMPLE  OUTPUT 

Table  3  presents  a  sample  output  listing  (TAPE6)  from  a  FREQRES 
calculation.  Record  number  4  (Figure  6)  was  the  input  data  record  and 
record  number  5  (Figure  7)  was  the  output  data  record.  The  output  echoes 
the  input  file  variable  specification.  This  output  listing  shows  that  the 
input  and  output  data  records  consisted  of  9970  points  each,  but  a  skip 
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factor  of  27  was  used.  The  large  skip  factor  was  used  in  order  to  force 
the  time  step  of  the  test  data  records  to  be  equal  to  the  time  step  of  the 
Speicher-Brode  waveform  created  from  a  FOURFIT  calculation.  The 
Speicher-Brode  waveform  time  step  (stated  three  fortfs  of  the  way  through 
the  output  listing)  is  1.337  x  10"^  seconds  and  the  original  test  data 
record  time  step  is  5  x  10“®  seconds.  The  skip  factor  of  27  now  gives 
the  test  data  records  a  time  step  of  (27) (5  x  10“^  seconds)  = 

1.35  X  10"^  seconds,  which  is  within  the  1  percent  tolerance  level  used 
in  the  program.  Next  the  listing  states  that  the  time,  input  data,  and 
output  data  conversion  factors  were  all  set  to  1.0.  The  beginning 
portions  of  the  input  and  Output  data  prior  to  2.67  msec  and  3.0  msec, 
respectively,  were  splined  to  zero.  The  final  15  percent  of  the  input  and 
output  data  were  also  splined  to  zero.  The  "very  important  notice"  in  the 
output  listing  compares  the  total  number  of  points  and  time  steps  between 
the  Speicher-Brode  waveform  and  the  test  data.  If  the  total  number  of 
points  is  not  identical,  then  the  program  truncates  the  time  history  with 
the  greatest  number  of  points  (the  Speicher-Brode  waveform  in  this 
example)  so  that  they  are  equal.  The  program  then  prints  out  the  time 
steps  for  user  inspection.  If  the  time  steps  do  not  meet  the  tolerance 
level  specified  in  FREQRES,  the  program  stops  and  the  final  three  lines 
shown  in  the  output  listing  will  not  be  printed.  If  the  tolerance  level 
is  met,  then  the  program  should  run  successfully  to  completion.  The  last 
three  lines  provide  the  value  of  the  full  integration  of  the  output  data 
record,  the  value  of  the  output  FFT  amplitude  spectrum  at  zero  frequency, 
and  the  value  of  the  offset  of  the  inverse  FFT  output  time  history  (see 
Section  2.3).  Notice  that  the  first  and  second  values  are  nearly 
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identical,  and  the  value  of  the  third  is  twenty  times  the  value  of  the 
second. 

Table  4  presents  the  input  file  variable  specifications  for  each  of 
the  test  data  records  as  used  in  this  analysis.  Only  two  of  the  data 
records  were  ever  used  as  input  data  records  in  the  cause-effect  analysis, 
record  numbers  2  and  4.  All  of  the  data  records  with  an  NSKIP  of  35  were 
used  as  output  data  records  when  record  number  2  v^as  used  as  the  input 
data  record.  All  of  the  data  recoros  with  an  NSKIP  of  27  were  associated 
with  record  number  4.  The  time  scale  was  always  kept  at  seconds  with 
TFAC  =  1.0.  Test  data  records  8  through  11  were  converted  from  units  of 
g's  to  ft/sec^  with  an  AFAC  =  32.2  and  then  integrated  so  that  the 
output  data  record  could  be  velocity  rather  than  acceleration  versus 
time.  All  of  the  strain  plots  were  inverted  with  a  AFAC  =  -1.0  to  be 
consistent  with  the  pressure  plots  in  which  compression  is  positive.  Test 
data  record  numbers  8  and  11  were  the  only  records  which  appeared  to  need 
baseline  correcting  as  noted  in  Table  4. 

Table  5  presents  the  total  impulse,  first  FFT  value,  and  inverse  FFT 
offset  for  test  data  records  4  through  20. 

B.4  PROGRAM  FREPLT 

Program  FREPLT  is  the  post  process  plotting  program  for  the  program 
FREQRES.  A  complete  listing  of  FREPLT  is  presented  in  Appendix  D.  FREPLT 
simply  reads  the  plot  file  created  from  a  FREQRES  run  (TAPE48)  and  creates 
eight  hard  copy  plots  on  the  DNA  Cyber  computer  at  Los  Alamos,  New 
Mexico.  The  eight  plots  are  as  follows: 

1.  Input  data  record  time  history 

2.  Output  data  record  time  history 

3.  Input  data  record  FFT  amplitude  spectrum  up  to  3000  Hz 
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4.  Output  cata  record  FFT  amplitude  spectrum  up  to  3000  Hz 

5.  FRF  amplitude  ratios  up  to  3000  Hz 

6.  Inverse  FFT  output  record  time  history 

7.  Output  record  FFT  with  new  input  data  record  influence 

8.  Modified  output  response  time  history 

B.5  PROGRAMMING  NOTES 

This  subsection  mainly  discusses  file  manipulation  among  the  three 
programs  FOURFIT,  FREQRES,  and  FREPLT.  First,  a  FOURFIT  calculation  is 
run  with  the  updates  discussed  in  Section  2.5  and  shown  in  Appendix  A. 
These  updates  create  TAPE49  which  must  be  saved. 

The  input  files  necessary  for  a  FREQRES  calculation  are  shown  in 
Table  6.  The  output  files  created  from  a  successful  FREQRES  calculation 
are  shown  in  Table  7.  FREPLT  only  requires  one  input  file,  TAPE9,  which 
is  the  plot  file  from  the  FREQRES  calculation,  TAPE48.  FREPLT  creates  two 
output  files,  TAPE6  and  PLOT.  TAPE6  is  printed  output.  PLOT  contains  the 
eight  plots  described  in  Secton  B  4  which  can  be  disposed  to  hard  copies 
with  a  job  control  instruction: 

PESP.0RIENT=R0TATE  MAJ0R="any  message" 

From  Table  6,  TAPE26  and  TAPE27  contain  digitized  test  data  records 
as  they  are  stored  on  EU  tapes  under  a  format  used  by  the  Waterways 
Experiment  Station  in  Vicksburg,  Mississippi.  The  format  is  readily 
apparent  in  the  listing  of  FREQRES  in  Appendix  C  with  the  READ  statements 
for  TAPE26  and  TAPE27. 

For  all  FFT's  and  FRF's,  the  maximum  frequency  of  concern  for 
plotting  is  assumed  to  be  3000  Hz.  There  is  very  little  significance  at 
frequencies  greater  than  3000  Hz  in  most  structural  and  soil  response 
characteristics.  Also,  the  power  in  most  FFT  amplitude  spectrums  of  test 
data  is  negligible  at  3000  Hz  and  beyond.  This  maximum  frequency  of 
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concern  (for  plotting  purposes  only)  can  be  altered  in  the  portion  of 
FREQRES  as  stated  below: 

C  Assumming  record  duration  is  .05  sec  and  maximum  frequency 
C  of  concern  is  3000  Hz  then  the  maximum  number  of  points  of  concern 
C  for  plotting  in  the  frequency  domain  is  {3000){.05)  =  150. 

NPFF  =  NPF 

IF{NPF.GT.150)  NPFF=150 

If  the  input  and  output  data  record  durations  are  different  from 
.05  seconds  or  if  the  maximum  frequency  of  concern  for  plotting  is 
different  from  3000  Hz,  then  the  150  in  the  last  statement  above  must  be 
altered  accordingly. 
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APPENDIX  C 


LISTING  OF  PROGRAM  FREQRES 


PROGRAM  FREQRES (INPUT, OUTPUT, TAPE5= INPUT, TAPE6-0UTPUT, 

*  TAPES,  TAPES, TAPES6, TAPES?, TAPE48) 

C 

It************************************** 

C 

C  PROGRAM  FREQRES  DETERMINES  THE  FREQUENCY  RESPONSE  FUNTION 
C  BETWEEN  TWO  SETS  OF  DATA  < "INPUT"  DATA  TRACE  X (NT) 

C  FROM  TAPE26  AND  "OUTPUT"  DATA  TRACE  A (NT)  FROM 

C  TAPES?).  THE  “ INPUT" /“OUTPUT"  LABELS  ARE  ONLY 
C  RELEVANT  IN  "CAUSE-EFFECT"  ANALYSIS  BETWEEN  TWO  SETS 

C  OF  DATA. 

C  RESULTS  ARE  WRITTEN  TO  A 

C  FILE  (TAPEAS)  TO  BE  READ  AND  PLOTTED  BY  PROGRAM 
C  FREPLT. 

C 

C*********************************************************************** 

C 

DIMENSION  FRQ (3000) , XFFT (3000) , AFFT (3000) ,  AMP (3000) , 

*TIMX (6000) , XTD (6000) , ATD (6000) , DDT ( 1S000) ,  PRESS (6000) 

DIMENSION  IX(8), IA(a),  ITX(6, 4), ITY(6, 4) 

DIMENSION  DUM(3), IWKE (5500) , WKE (5500) 

EQUIVALENCE  ( IWKE ( 1 ), WKE ( 1 ) > 

DIMENSION  COPS (3000), Cl PS (3000) 

COMPLEX  XFFT, AFFT 
C 

c 

C  TAPES  INPUT  FILE  DESCRIPTION 

C 


c 

p 

CARD 

COLUMN 

FORMAT 

VARIABLE  DESCRIPTION 

w 

c 

1 

1-5 

15 

NEPTS 

NO.  OF  POINTS  TO  BE  READ  FROM  TAPE 

c 

6-10 

15 

NSKIP 

SKIP  INTERVAL  (DEFAULT=1,  ALL  POINTS 

c 

FROM  TAPE  ARE  SAVED) 

c 

11-20 

E10.  3 

TFAC 

TIME  CONVERSION  FACTOR  (DEFAULT=1. 0) 

c 

21-30 

E10.  3 

XFAC 

INPUT  DATA  CONV.  FACTOR  (DEFAULT=1 . 0) 

c 

31-40 

E10.  3 

AFAC 

OUTPUT  DATA  CONV.  FACTOR  (DEFAULT=1 . 0) 

C 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


i-i 


6-lS 


16-20 


£1-30 


31-35 


36-45 


46-50 


51-60 


1-5 


6-15 


15 


ISPBX 


E10. 3  TSPBX 


15 


ISPEX 


E10. 3  TSPEX 


15 


ISPBA 


El  0.3  TSPBft 


15 


ISPEPI 


El  0.3  TSPEfi 


15 


IBLX 


El  0.3  DELPX 


0: 


1: 


NO  SPLINE  PERFORMED  ON  BESINNING 
OF  INPUT  DflTft 

BEGINNING  OF  INPUT  DfiTft  WILL  BE 
SPLINED  (DEFflULT=0) 

IF  ISPBX=1,  A  SPLINE  IS  PERFORMED  AT 
THIS  TIME  BACK  TO  TIME  ZERO  <DEFAULT 
IS  0.0  WHICH  MEANS  NO  SPLINE  IS  DONE) 
0:  NO  SPLINE  PERFORMED  ON  END  OF 
INPUT  DATA 

1:  END  OF  INPUT  DATA  WILL  BE  SPLINED 
(DEFAULT=0) 

IF  ISPEX=1,  TIME  AT  WHICH  SPLINE 
BEGINS  FOR  INPUT  DATA  (DEFAULT  IS 
8554  OF  TTOT) 

0:  NO  SPLINE  PERFORMED  ON  BEGINNING 
OF  OUTPUT  DATA 

1:  BEGINNING  OF  OUTPUT  DATA  WILL  BE 
SPLINED  (DEFAULT=0) 

IF  ISPBA=1,  A  SPLINE  IS  PERFORMED  AT 
THIS  TIME  BACK  TO  TIME  ZERO  (DEFAULT 
IS  0.0  WHICH  MEANS  NO  SPLINE  IS  DONE) 
NO  SPLINE  PERFORMED  ON  END  OF 
OUTPUT  DATA 

END  OF  OUTPUT  DATA  WILL  BE 
SPLINED  (DEFAULT=0> 

ISPEA=1,  TIME  AT  WHICH  SPLINE 
BEGINS  FOR  OUTPUT  DATA  (DEFAULT 
IS  8554  OF  TTOT) 


0: 


1: 


IF 


INPUT  DATA  BASELINE  CORRECTION  TRIGGER 
0:  NO  BASELINE  CORRL;TION 
1:  BASELINE  CORRECTION  WITH  THE 

FOLLOWING  PARAMETERS  (DEFAULT=0) 
CORRECTION  ADDED  TO  INPUT  DATA  VALUES 
AFTER  TIME  SBX.  IF  EBX  AND  SBX  ARE 
EQUAL  THEN  THE  FULL  VALUE  OF  DELPX  IS 
ADDED  AT  ALL  TIMES  AFTER  SBX.  IF  EBX 
IS  GREATER  THAN  SBX  THEN  THE  PLOT  IS 
ROTATED  ABOUT  THE  POINT  DEFINED  AT  SBX 
BY  THE  AMOUNT  DELPX  AT  TIME  EBX. 


c 

c 

16-25 

E10.  3 

SBX 

START  TIME  FOR  BASELINE  CORRECTION 
(PLOT  ROTATION  POINT  IF  EBX  >  SBX) 

c 

£6-35 

E10.  3 

EBX 

END  TIME  FOR  BASELINE  CORRECTION 

c 

c 

c 

36-40 

15 

INTX 

INPUT  DATA  INTEGRATION  TRIGGER 

0:  NO  INTEGRATION 

1;  INTEGRATE  INPUT  DATA 

c 

C  4 

C 

1-5 

15 

IBLA 

OUTPUT  DATA  BASELINE  CORRECTION  TRIGGER 
0:  NO  BASELINE  CORRECTION 

1:  BASELINE  CORRECTION  WITH  THE 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


£-15 


E10. 3  DELPft 


NOTE: 


FOLLOWING  PPlRftMETERS  <DEFPULT=0) 
CORRECTION  ODDED  TO  OUTPUT  DftTfi  VPLUES 
AFTER  TIME  SBA.  IF  EBA  AND  SBA  ARE 
EQUAL  THEN  THE  FULL  VALUE  OF  DELPA  IS 
ADDED  AT  ALL  TIMES  AFTER  SBA.  IF  EBA 
IS  GREATER  THAN  SBA  THEN  THE  PLOT  IS 
ROTATED  ABOUT  THE  POINT  DEFINED  AT  SBA 
BY  THE  AMOUNT  DELPA  AT  TIME  EBA. 

START  TIME  FOR  BASELINE  CORRECTION 
(PLOT  ROTATION  POINT  IF  EBA  >  SBA) 

END  TIME  FOR  BASELINE  CORRECTION 
OUTPUT  DATA  INTEGRATION  TRIGGER 
0:  NO  INTEGRATION 
1:  INTEGRATE  OUTPUT  DATA 


ALL  OF  THE  FOLLOWING  LABELS  SHOULD  BE  CENTERED  WITHIN  THE 
FIRST  30  COLUMNS  OF  EACH  LINE  OF  THE  INPUT  FILE. 


16-^5 

E10.  3 

SBA 

£6-35 

E10.  3 

EBA 

36-40 

15 

I  NT  A 

5 

1-40 

4A10 

ITX 

INPUT  DATA  X-AXIS  LABEL;  EXAMPLE: 
TIME  (SEC) 

6 

1-40 

4A10 

ITY 

INPUT  DATA  Y-AXIS  LABEL;  EX: 
PRESSURE  (PSD 

7 

1-40 

4A10 

ITX 

OUTPUT  DATA  X-AXIS  LABEL;  EX: 

TIME  (SEC) 

8 

1-40 

4A10 

ITY 

OUTPUT  DATA  Y-AXIS  LABEL;  EX: 
STRAIN  (IN/IN) 

9 

1-40 

4A10 

ITX 

INPUT  DATA  FOURIER  AMP.  SPEC. 
X-AXIS  LABEL;  EX: 

FREQUENCY  (HZ) 

10 

1-40 

4A10 

ITY 

INPUT  DATA  FOURIER  AMP.  SPEC. 
Y-AXIS  LABEL;  EX: 

AMPLITUDE  (PS I -SEC) 

1 1 

1-40 

AA10 

ITX 

OUTPUT  DATA  FOURIER  AMP.  SPEC. 
X-AXIS  LABEL;  EX: 

FREQUENCY  (HZ) 

1£ 

1-40 

4A10 

ITY 

OUTPUT  DATA  FOURIER  AMP.  SPEC. 
Y-AXIS  LABEL;  EX: 

AMPLITUDE  (SEC) 

13 

1-40 

4010 

ITX 

FREQUENCY  RESPONSE  FUNCTION 

X-AXIS  LABEL;  EX: 

FREQUENCY  (HZ) 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


14 


15 


16 


1-40  4fil0  ITY 


1-40  4fti0  ITX 


FREQUENCY  RESPONSE  FUNCTION 
Y-ftXIS  LABEL;  EX: 

FFTA/FFTX 

ERODE  OUTPUT  RESPONSE 
X-AXIS  LABEL;  EX: 

TIME  (SEC) 


1-40  4A10  ITY  BRODE  OUTPUT  RESPONSE 

Y-AXIS  LABEL;  EX* 
STRAIN  (IN/IN) 


C****-*********-****************^*^***************************************** 


c 

C  IF  YOU  ENCOUNTER  A  CM  LIMIT  ERROR  THEN  MOST  LIKELY  THE  ARRAY  SIZES 
C  FOR  WKE  AND  IWKE  SPECIFIED  ARE  TOO  SMALL  FOR  THE  VALUE  OF  NPT 

C  PASSED  THROUGH  THE  CALL  TO  FFTRC  OF  FFTCC.  THIS  ERROR  CAN  BE 

C  AVOIDED  BY  INCREASING  THE  SIZE  OF  THE  WKE  AND  IWKE  ARRAYS  AND  AT 

C  THE  SAME  TIME  INCREASING  THE  CM=?  AND  RFL, ?  SPECIFICATIONS  IN  THE 

C  JCL.  NOTE:  FOR  A  5500  SIZE  OF  WKE  AND  IWKE,  ?=S70000. 

C 

C************************************  ************************************* 

c 

REWINDS 
REWINDS 
WRITE (6, 199) 

199  FORMAT (*  FREQRES  OUTPUT  LISTING*) 

READ  (£,111)  NEPTS,  NSKIP,  TFAC, XFAC, AFAC 
111  rORMAT(£I5, 3E10. 3) 

SKIP=FLOAT(NSKIP) 

IF  (TFAC.  EQ.  0.  )  TFAC=1.0 
IF  (XFAC.  EQ.  0.  )  XFAC=1.0 
IF  (AFAC.  EQ.  0.  )  AFAC=1.0 

READ  (£,  1 15)  ISPBX, TSPBX,  ISPEX,  TSPEX,  ISPBA, TSPBA, ISPEA, TSPEA 
115  F0RMAT(4(I5,  E10.  3) ) 

READ(£, 1 1£)  IBLX,  DELPX,  SBX, EBX, INTX 
READ(£, 1 1£)  IBLA, DELPA,  SBA, EBA, INTA 
11£  FORMATdS,  3E10.  3,  15) 

113  FORMAT (4A10) 

DO  10  1=1,6 

READ(£,113)  (ITXd,  J),  J=l,4) 

10  READ(£,  113)  (ITYd,  J),  J=l,4) 

11  CONTINUE 

WRITE(6, 1£7)  NEPTS, NSKIP 
WRITE (6, 123)  TFAC 
WRITE(6, i£9)  XFAC 
WRITE (6, 130)  AFAC 

IF  (ISPBX. EQ.l)  WRITE (6, 131)  TSPBX 

IF  dSPEX.  EQ.  1.  AND.  TSPEX.  LE.  0.  0)  WRITE  (6,  132) 
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IF  (ISPEX. EQ. 1. ftND. TSP£X.  BT.  0.  0)  WRITE (6, 1/3)  TSPEX 
IF  <ISPBA. EQ-1)  WRITE <£, 134)  TSPBP 
IF  (ISPEft.EG. l.flND.TSPEft.  LE. 0. 0)  WRITECS, 135) 

IF  <ISPEA. EQ. 1. fiWD.  TSPEft.  GT.  0.  0)  WRITE'S, 136)  TSPEfl 
IF  (2BLX.EQ.1)  WRITE  <S,1£3)  DELPX,  SSX,  EBX 
IF  (IBLft. EQ. 1)  WRITECS, ISA)  DELPft, B3A, EBP 
IBS  FORMATC/,*  BASELINE  CORRECTICN  REQUESTED  FOR  INPUT  DATA  RECORD* 
+,*  WITH:*,/,*  DELPX  =  *,£10.3,*;  SBX  =  *,E10.3,*;  EBX  =  * 
+,E10.  3) 

IBA  FORMATC/,*  BASELINE  CORRECTION  REQUESTED  FOR  DLTPLT  DATA  RECORD* 
+,*  WITH;*,/,*  DELPA  =  *,E10.  3,  *;  SBA  =  *,E10.3,  *;  E3A  =  * 
+,E10.  3) 

137  FORMATC/,*  THE  NUMBER  OF  POINTS  READ  FROM  THE  DATA  RECORD  TAPES* 
+,*  IS  *,  15,  *  WITH  A  SKIP  OF  *,!£,*  CONSIDERED  FOR  ANALYSIS*) 

128  FORMAT </,*  TIME  CONVERSION  FACTOR  =  *,E10.  3) 

129  FORMAT <*  INPUT  DATA  CONV.  FACTOR  =  *,  E10.  3) 

130  FORMAT (*  OUTPUT  DATA  CONV.  FACTOR  =  *,  E10.  3) 

131  FORMAT (/,*  BEGINNING  PORTION  OF  INPUT  DATA  SPLINED  FROM  TIME* 

+,*  EQUAL  0.0  TO  TIME  EQUAL  *,E10.  3) 

132  FORMAT  (/,*  FINAL  15)4  OF  INPUT  DATA  SPLINED  TO  ZERO*) 

133  FORMAT </,*  INPUT  DATA  FROM  TIME  =  *,E10. 3, *  ON  SPLINED  TO  ZERO*) 

134  FORMAT </,*  BEGINNING  PORTION  OF  OUTPUT  DATA  SPLINED  FROM  TIME* 
+,*  EQUAL  0.0  TO  TIME  EQUAL  *,  ElO.  3) 

135  FORMATC/,*  FINAL  IS'/.  OF  OUTPUT  DATA  SPLINED  TO  ZERO*) 

136  FORMATC/,*  OUTPUT  DATA  FROM  TIME  =  *, ElO. 3, *  ON  SPLINED  TO  ZERO*) 
WRITE  <48, 111)  NEPTS,  NSKIP,  TFAC, XFAC, AFAC 

WRITE  <48,  1 IS)  ISPBX, TSPBX,  ISPEX, TSPEX, ISPBA,  TSPBA,  ISPEA, TSPEA 
WRITE <48, 112)  IBLX,  DELPX, SBX, EBX 
WRITE  <48, 1 12)  IBLA, DELPA, SBA, EBA 
REWIND2S 
REWIND27 
C 

C  READ  INPUT  DATA  RECORD  TAPE  HEADER  INFORMATION 
C 

READ<£6,  30)  IX<3),IX<4), 


* 

DUM<1),DUM<2) 

* 

IX<1),  IX<2), 

* 

DTD, NP 

30  FORMAT<3<2A10) , ElS.a, 15) 

IF  <E0F<26))  900,901 
901  NPT=NEPTS/NSKIP 
DTD=TFAC*SK I P*DTD 
TIMX<1)=0. 

DO  12  1=2, MPT 
12  TIMX <I)=TIMX <I-1)+DTD 
TTOT=TIMX<NPT) 

IF  <TSPBX  LT. 0. 001*TTOT. OR. TSPBX. GT. TTOT)  ISPBX=0 
IF  <TSPBA. LT. 0. 001*TTOT. OR. TSPBA. GT. TTOT)  ISPBA=0 
C 

C  READ  INPUT  DATA  RECORD  VALUES 
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READ (£6, 50)  (DDT < I) , 1=1 , NEPTS) 

50  FORMAT (5E IS. 8) 

IF  (EOF (£6))  300, 90£ 

90E  IM=0 

DO  13  I=NSKIP, NEPTS, NSKIP 

IM=IM+1 

XTD<IM)=0. 

DO  14  J=l,NSKiP 

14  XTD( iM)=XTD(IM)+DDT(I+l-J) 

13  X‘1D(IM)  =  (XTD(IM)  /SKIP)  *XFAC 

IP"  <IM.  LT.NPT)  \’PT=IM 
C 

C  BASELINE  COR=?£CTION  FOR  INPUT  DATA 
C 

IF  (IBLX.EQ.  0)  GO  TO  £1 
DO  15  1=1, NP7 

15  IF  (SBX. LT. TIMX (I ) ^  GO  TO  18 

16  ISBX=I 

IF  ( (EBX-SBX) . GT. (0. 00i*TTOT) )  GOTO  101 
IF  (ISBX. GE. NPT)  GO  TO  9S0 
GO  TO  10£ 

101  CONTINUE 

DO  17  I  =  ISBX,Ni='T 

17  IF  (EBX.lT.TIMXO)  go  TO  18 

18  IEBX=I 

IF  (lEBX. LE. ISBX. OR. lEBX. GE. NPT)  GO  TO  990 
DO  19  I=ISBX,NPT 

19  XTD(I)=XTD(I)+DELPX* ( ( I-ISBX) / ( lEBX-ISBX) ) 

GO  TO  £1 

10£  DO  £0  I=ISBX,NPT 
£0  XTD(I)=XTD<I)+D£LPX 

21  IF  (INTX.EQ. 1)  CALL  I NT6RT (DTD, NPT, TIMX, XTD, DDT) 

IF  (ISPEX.EQ.  1)  CALL  SPLINE (TS3EX, TTOT,  NPT,  TIMX, XTD) 
IF  (ISPBX.EQ.l)  CALL  BS^'LIN  (TSPBX,  TTO^,  N^'T,  TIMX,  XTD) 
IX<5)=10H  INPUT  DAT 
IX(6)=10HA  RECORD; 

IX (7)=10HX (NT) 

IX (3)=10H 

call  PLTSAVd,  IX,  ITX,  ITY,  NPT,  TIMX,  XTD) 

C 

C  READ  OUTPUT  DATA  RECORD  TAPE  HEADER  INFORMATION 
C 

READ<£7, 30)  IA(3),IA(4), 

*  DUM(l) , DUM(£) , 

*  IA(1) , IA(£) , 

*  DAM, NP 
IF  (E0F^£7))  900,903 

903  CONTINUE 
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C  READ  OUTPUT  DATA  RECORD  VALUES 
C 

READ <£7, 50)  (DDT < I > ,  1  =  1,  NEPTS) 

IF  (EOF (£7))  900,904 
904  IM=0 

DO  £3  I=NSKIP,  hJEPTS,  NSKIP 

IM=IM+1 

ATD(IM)=0. 

DO  £4  J=1,NSKIP 
£4  ATD(IM)=ATD(IM)+DDT(I+3-J) 

£3  ATD(IM)=(ATD(IM) /SKIP)*AFAC 
C 

C  BASELINE  CORRECTION  FOR  OUTPUT  DATA 
C 

IF  (IBLA. EQ. 0)  GO  TO  31 
DO  £5  1=1, NPT 

£5  IF  (SBA. LT. TIMX (I ) )  GO  TO  £6 
£6  ISBA=I 

IF  ( (EBA-SBA) . GT.  (0.  001*TTGT) )  GO  TC  103 
IF  (ISBA. GE. NPT)  GO  TO  991 
GO  TO  104 

103  CONTINUE 

DO  £7  I=ISBA,NPT 
£7  IF  (EBA. LT. TIMX (I) )  GO  TO  £8 
£8  IEBA=I 

IF  (lEBA.  LE.  ISBA.  OR.  lEBA.  GE.  NPT)  GOTO  991 
DO  £9  I=ISBA,NPT 

29  ATD ( I ) =ATD ( I ) +DELPA* ( ( I-IS3A) / ( lEBA-ISBA) ) 

GO  TO  31 

104  DO  39  I=ISBA,NPT 

39  ATD ( I ) =ATD ( I ) +DELPA 

31  IF  (INTA.EQ. 1)  CALL  INTGRT (DTD, NPT, TIMX, ATD, DDT) 

IF  (ISPEA.EQ. 1)  CALL  SPLINE (TSPEA, TTOT, NPT, TIMX, ATD) 

IF  (ISPBA.EQ.  1)  CALL  BSPLIN(TSPBA, 'T'TOT,  NPT,  TIMX,  ATD) 
IA(5)=10H  OUTPUT  DA 
IA(6)=10HTA  RECORD: 

IA(7)=10H  A(NT) 

IA(a)=10H 

CALL  PLTSAV (£, lA, ITX, ITY, NPT, TIMX, ATD) 

SUMIMP=0. 

DO  80  1=2, NPT 

AREA= (ATD ( I-l ) +ATD ( I ) ) *DTD/£. 

80  SUMIMP=SUMIMP+AREA 
C 

C  READ  THE  SPEICHER-BRODE  FIT  TO  THE  INPUT  DATA  FROM  TAPE3 

C 

READ(3,£00)  NBPTS,DTBP 

READ(3,£10)  TIMl,  ( (TIMD),  IU=£, NBPTS) 

READ (3, £10)  (PRESS ( IP) , IP=1, NBPTS) 

200  r  URMATdS,  E15.  8) 
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nnnnn  non 


£10  FORMfiT(10E15.a) 

IF  (TIMl.LE.  0.  >  GO  TO  £51 
NPLUS=TIM1/DTEP 
NBPTS=NBPTS+NPLUS 
NEND=NBPTS~NPLUS 
DO  £5£  IP=1,NEND 

£5£  PRESS  <NBPTS+1-IP) =PRESS  <NBPTS-NPLUS-IP+1 ) 

DO  £53  IP=1,NPLUS 
£53  PRESS ( IP) =0. 

C  CONVERT  THE  ERODE  FIT  FROM  MPA  TO  PS I 
£51  DO  £54  IP=1,NBPTS 
£54  PRESS (IP) =PRESS ( IP) *145. 0377 
WRITE <&,  441) 

441  FORMAT(//, 1X,80<1H*)//) 

WRITE (&,44£)  NBPTS, NPT 
44£  FORMAT (*  VERY  IMPORTANT  NOTICE:*/ 

**  NBPTS  FOR  THE  SPEICHER-BRODE  FIT  FROM  FDORFIT  =  *, 15/ 

**  MUST  BE  EQUAL  TO  NPT  FROM  THIS  PROGRAM  =  *,15/ 

**  IF  NBPTS  AND  NPT  ARE  NOT  EQUAL,  THIS  PROGRAM  WILL  */ 

**  TRUNCATE  ONE  OF  THEM  TO  MAKE  THEM  EQUA_. */) 

WRITE <6, 443)  DTBP, DTD 

443  FORMAT(*  DTBP  FOR  THE  SPEICHER-BRODE  FIT  FROM  FOURFIT  =*,E15. 8/ 

*  *  MUST  BE  VERY  CLOSE  TO  DTD  FROM  THIS  PROGRAM  =*, E15. 8/ 

**  IF  THEY  ARE  NOT  CLOSE,  THIS  PROGRAM  WILL  STOP.*/) 

WRITE <6, 441) 

IF  (NPT. GT. NBPTS)  NPT=NBPTS 

IF  (DTBP.  LT.0.99*DTD.  OR.  DTBP.  GT.  1.01*DTD)  STOP 

FOURIER  TRANSFORM  OF  INPUT  DATA  RECORD 

DFE=1. /TTOT 
FRa(l)=0. 0 

CALL  FFTRC ( XTD, NPT,  XFFT,  IWKE, WKE) 

NPF«NPT/£+1 

ASSUMMING  RECORD  DURATION  IS  .05  SEC  AND  MAXIMUM  FREQUENCY 
OF  CONCERN  IS  3000  HZ  THEN  THE  MAXIMUM  «  OF  POINTS  OF  CONCERN 
FOR  PLOTTING  IN  THE  FREQUENCY  DOMAIN  IS  (3000) (.05)  =  150. 

npff=npf 

IF  (NPF.GT. 150)  NPFF=150 
XRE=REAL ( XFFT ( 1 ) ) /NPT 
XIE=AIMAG(XFFT(1) ) /NPT 
AMP ( 1 ) =SQRT ( X  RE**£+X I E**£ ) *TTOT 
DO  3£  JK=£,NPF 
FRQ ( JK) =FRQ ( JK-1 ) +DFE 
XRE=REAL(XFFT(JK) ) /NPT 
XIE=AIMAG(XFFT(JK) ) /NPT 
3£  AMP ( JK ) =SQRT ( X RE**£+X I E**£) *TTOT 
IX(5)=10H  FOURIER  A 
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nnn  on  nnnn  cinn 


IX (6)=10HMPLITUDE  S 
IX  (7)=^10HPECT1RUM;  X 
IX  (a)=10H(^4/TTOT) 

COLL  PLTSflV(3,  IX,  ITX,  ITY,  NPFF,  FRQ,  PiMP) 

FOURIER  TRANSFORM  OF  OUTPUT  DftTft  RECORD 

CALL  FFTRC  (PTD,  NPT,  PiFFT,  IWKE,  WKE) 

DO  33  JK=1,NPF 
ARE=REAL (AFFT ( JK) ) /NPT 
AIE=fti:^fiG<AF.-T(  JK)  )  /NPT 

33  AMP ( JK) =SQRT ( ARE**£+AIE**£) *TTOT 
FFT1=AMP(1) 

IA<5)=10H  FOURIER  A 
IA<6)=12iHMP_ITUDE  S 
IA<7)=10:"iPEOTRUM:  A 
IA(a)=10l-!<.\/':TCT) 

CALL  PLTGAV  <4,  lA,  ITX, ITY, NPFF, FRQ, AMP) 

CALCULATE  THE  FREQUENCY  RESPONSE  FUNCTION  (COMPLEX) 

AND  PLOT  UP  ITS  MAGNITUDE  (AMPLITUDE  RESPONSE  RATIO) 

DO  34  JK=1,NPF 
XRE=REAL(XFFT(JK) ) /NPT 
XIE=AIMAG(XFFT(JK) ) /NPT 
ARE=REAL (AFFT ( JK) ) /NPT 
AIE=AIMAG(AFFT(JK) > /NPT 

COPS(JK)=(ARE*XRE+AIE*XIE) / (XRE**£+XIE**£) 

COPS(JK)  IS  THE  REAL  PART  OF  THE  FREQUENCY  RESC'ONSE  FUNCTION 
CIPS(  JK)  =  (XRE*AIE-ARE*XIE)  /  ( XRE«-*£+X IE**£) 

CIPS(JK)  IS  THE  IMAGINARY  PART  OF  T.-iE  FREQUENCY  RESPONSE  FUNCTION 
AMP  (JK)==SQRT  (COPS  (JK)  **£+CIPS(  JK)  **£) 

34  CONTINUE 
IX (1)=10H 

IX(£)=10H  FRE 

IX (3)=10HQUENCY  RES 
IX (4)=10HPONSE  FUNC 
IX (5)=10HTION  (AMPL 
IX (6)=10HITUDE  RATI 
IX (7)=10HOS) 

IX (a)=10H 

CALL  PLTSAV(5, IX, ITX, ITY, NPFF, FRQ, AMP) 

INVERSE  FFT  FOR  OUTPUT  DATA  RECORD 

DO  81  1=1, NPF 

AFFT ( I ) =CON JG ( AFFT ( I ) ) 

81  CONTINUE 

CALL  FFTCC ( AFFT, NPF, IWKE, WKE) 

DO  8£  1=1, NPF 
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«FFT  < I ) =CDN JG  <  flFFT ( I ) ) /NPF 
8£  CONTINUE 

DTD^T I MX ( NPT ) / ( NPF- 1 ) 

DO  83  1=1, NPF 
TIMX(I+1)=TIMX (I)+DTD 
ftRE=REflL(PFFT (I> ) 
flIE=ftIMflG(ftFFT(I) > 

83  flTD(I)=fiRE 
OFFSET=PlTD  (1) 

DO  84  1=1, NPF 

84  AID ( I ) =PTD ( I ) -OFFSET 

IA(5)=:2H  OUTPUT  RE 
Ifi  (&)=iZ:-;CC:RD  FROM 
Ift  (7)  =15?F INVERSE  FF 
lA  <8)  " 

CALL  PLTSAV (£,  IR,  ITX,  ITY, NPF, TIMX, AID) 

WRITE <6, 444)  SUM IMP, FFTl,  OFFSET 
444  raRMAT(*  OUTPUT  RECORD  TOTAL  IMPULSE  =  *,  E15.&/ 

*  FIRST  POINT  OF  OUTPUT  FFT  =  *,  E15.  6/ 

*  *  OUTPUT  INVERSE  FFT  OFFSET  =  *,E15. &//) 

C 

C  FOURIER  TRANSFORM  OF  THE  SPEICHER-BRODE  FIT 
C 

CALL  FFTRC (PRESS,  NPT,  XFFT,  IWKE,  WKE) 

DO  59  JK=1,NPF 
XRE=REAL(XFFT(JK) )  'NPT 
XIE=f>I.MAB(XFFT(JK)  )  /NPT 
COPSO--COPS  ( JK )  *X  RE-C I  PS  ( JK )  *  X I E 
C I PSD=COPS ( JK ) *X I E+C I PS ( JX) *X  RE 
COPS ( JK) =COPSD*NPT 
59  CIPS<JK)=CIPSD*NPT 

C  COPS(JK)  IS  NOW  THE  REAL  PART  OF  THE  FFT  OF  THE  OUTPUT  RESPDNS 

C  IF  THE  SPEICHER-BRODE  FIT  WOULD  HAVE  BEEN  -^HE  LOADING  INPUT 

C  CIPS(JK)  IS  NOW  THE  IM^/GINARY  PART  OF  THE  F-T 
DO  40  JK=r,NPF 

40  AFFT<JK)=CMPLX (COPL  i  () , CIPS ( JK) ) 

DO  45  JK=1,NPF 
ARE=REAL (AFFT ( JK) ) /NPT 
AIE=AIMRG(AFFT(JK) ) /NPT 
45  AMP  ( JK)  =SQRT  ( ARE**£'+AIE**£)  *TTOT 
IA(5)=10HSP-BRODE  F 
IA(6)=10HFT  AMPLITU 
IA(7)=10HDE  SPECTRU 
IA(8)=10HM 

CAL.  PLTSAV (4,  I A,  ITX,  ITY,  NPFF,  FRQ, AMP) 

C 

C  INVERSE  FFT  THE  SPEICHER-BRODE  OUTPUT  RESPONSE  RFFT(JK) 

C 

DO  41  1=1, NPF 

AFFT  < I ) =CON JG ( AFFT ( I ) ) 
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c 

c 

c 

c 


41  CONTINUE 

CALL  FFTCC (AFFT, NPF, IWKE, WKE) 
DO  42  1=1, NPF 

AFFT  < I ) =CONJG (AFFT  < I ) ) /NPF 

42  CONTINUE 
ARE=REAL<AFFT(I) ) 
AIE=AINAG(AFFT<I) ) 

43  ATD(I)=ARE 
OFFSET=ATD ( 1 ) 

DO  44  1=1, NPF 

44  ATD(I)=ATD(I)-OFFSET 
IA<5)=10H  OUTPUT  RE 
IA(&)=liZiHCORD  IF  LO 
IA(7)=12iHADING  WERE 
IA(a)=l!ZlH  SP-BRODE 


CALL  PLTSAV (S, lA, ITX, ITY, NPF, TIMX, ATD) 

GO  TO  999 
900  WRITE (6, 70) 

70  FORMAT  (i:-l  ,  !^SND-0F-F1LE  REACHED  EARLY-^, ///) 

GO  TO  999 

990  WRITE(£,71) 

71  FORMAT (IH  , *ERROR  IN  BASELINE  CORRECTION  INSTRUCTIONS  FOR 
+*INPUT  DATA*,///) 

GO  TO  999 

991  WRITE <6, 72) 

72  FORMAT <1H  , *ERROR  IN  BASELINE  CORRECTION  INSTRUCTIONS  FOR 
+*OUTPUT  DATA*,///) 

999  END 

SUBROUTINE  PLTSAV (N,  ITLT,  I  XL, lYL, NN, XP, YP) 

DIMENSION  ITLTd),  I  XL  <6,  1),  lYLCG,  1 ) ,  XP  ( 1 ) ,  YP  ( 1 ) 

WRITE (48, 35)  ( ITLT (L) ,  L=l, 8) 

35  FORMAT (aA10) 

WRITE(4e,36)  (IXL(N, L) , L=l, 4) 

36  FORMAT (4A10) 

WRITE(48,36)  ( lYL (N, L) , L=l, 4) 

CALL  FMAX (XP, NN. XPMN, XPMX) 

CALL  FMAX (YP, NN,  YPMN, YPMX) 


WRITE (48, 37)  NN, XPMN, XPMX, YPMN, YPMX 

37  F0RMAT(I5,4E15. 8) 

WRITE  (48, 38)  ( XP (K) , K=1 , NN) 

WRITE (48, 38)  (YP (K) , K=l, NN) 

38  FORMAT (10E15. 8) 

RETURN 

END 

SUBROUTINE  FMAX (ARY, NA, XMN, XMX ) 


*******1h*¥-************************  ******************** 


THIS  SUBROUTINE  FINDS  THE  MAXIMUMS  AND  MINIMUMS 
OF  THE  VARIOUS  ARRAYS  TO  BE  PLOTTED  BY  FREPLT 

***^^**•^^•^^*■)^■■){•**•^H^•^^**«•*^<■****^<•*^(•*****■)(•***^(•*********■^<•*i<■*** 

DIMENSION  ARY(l) 
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c 

XMN  =  ftRY(l) 

XMX  =  ftRY<l) 

IF(NA. EQ. 1)  RETURN 
DO  10  1=2, NA 

IF(XMN.  GT.  ARY(I)  )  XNN  =  ARYC) 

IF(XMX.LT.  ARY<I)  )  XNX  =  ARY<I) 

10  CONTINUE 
C 

RETURN 

END 

SUBROUTINE  SPLINE (TBGN,  TLAST,  NP, TTIM, PRESS) 
C********************************************************** 
C  THIS  SUBROUTINE  SETS  UP  A  COSINE  SOUARED  SPLINE 
C  FUNCTION  AND  APPLIES  IT  TO  THE  FINAL  PORTION 
C  (DEFAULT=15-/4)  OF  THE  DATA  RECORD  TO  AVOID  AN 
C  INFINITE  IMPULSE  AT  ZERO  FREQUENCY  FOR  RECORDS 
C  WHICH  DO  NOT  RETURN  TO  ZERO 

C******************************** '<■************************* 
DIMENSION  TTIM<1),PRESS(1) 

C 

PIE=3. 1415927 
SPN=TBGN/TLAST 
IF  (SPN.LT.  0.01)  SPN=0.85 
K=IFIX<SPN*NP) 

IM=IMP-K+1 

T1=TTIM<K) 

DO  10  J=1,N 

TFACT=  <TTIM (K) -T1 ) / <TLAST-T1 ) 

SFACT=COS (TFACT*PIE*.  5) 

SFACT=SFACT**2 

PRESS (K) =PRESS (K) *SFACT 

K=K+1 

10  CONTINUE 
RETURN 
END 

SUBROUTINE  BSPLIN(TBGN,  TLAST, NP,  TTIM,  PRESS) 
C******************************************************** 

C  THIS  SUBROUTINE  APPLIES  A  COSINE  SQUARED  SPLINE 
C  TO  THE  BEGINNING  PORTION  OF  THE  DATA  RECORD  TO 
C  AVOID  AN  INFINITE  IMPULSE  AT  ZERO  FREQUENCY  FOR 
C  RECORDS  WHICH  DO  NOT  START  AT  ZERO 

C************************************-*^t****************** 

c 

DIMENSION  TTIM(1),PRESS(1) 

C 

PIE=3. 1415927 
SPN=TBGN/TLAST 
K=IFIX (SPN*NP) 

T1=TTIM(K) 
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DO  10  J=1,K 
TFftCT=<Tl-TTlM<J> )/Tl 
SFOCT«COS  <TFftCT*PIE*. 5) 

SFPlCT=SFftCT**2 
10  PRES5  (j) —PRESS  (J)  ■kSFACT 
RETURN 
END 

SUBROUTINE  INTGRT  <DTD, NP, TTIM, PRESS, DUM) 

C  THIS  SUBROUTINE  PERFORMS  AN  INTEGRATION  OF 
C  DIGITIZED  DATA  ACCORDING  TO  THE  TRAPEZOIDAL  RULE 
C******************************************************* 
DIMENSION  TTIM  < 1 > , PRESS  < 1 ) , DUM  < 1 ) 

PRESS  < 1 ) »0. 

SUMIMP-0. 

DO  80  J=£, NP 

AREA" <PRESS  < J-1 ) +PRESS ( J) ) *DTD/£. 

SUM I MP"SUM I MP+ AREA 

80  DUM<J)"SUMIMP 
DO  81  I=£,NP 

81  PRESS<I)=DUM(I> 

RETURN 

END 


APPENDIX  D 


LISTING  OF  PROGRAM  FREPLT 


PROGRftM  FREPLT (INPUT,  OUTPUT,  TftPE5=INPUT, TflPE6=0UTPUT, 

*  TAPES, PLOT) 

COMMON  /PLOTV/  ITL (8) , ISTL (8) , IDB 
COMMON  /PID/  NPTS,NSKIP,TFAC, XFAC, AFAC, 

*ISPBX, TSPEX, ISPEX,  TSPEX,  ISPBA, TSPBA, ISPEA, TSPEA, 

*IBLX, DELPX, SBX, EBX,  IBLA,  DELPA, SBA, EBA 
DIMENSION  XARY (9000) ,  VARY (9000) , LABX (A) , LABY (4) 

CALL  GPLOTdHU,  7HARABINB,  7) 

CALL  BGNPL<-1) 

READ (9, 111)  NPTS,NSKIP, TFAC, XFAC, AFAC 

READ (9, 115)  ISPBX, TSPEX, ISPEX, TSPEX, ISPBA, TSPBA, ISPEA, TSPEA 
READ (9, 112)  IBLX,  DELPX, SBX, EBX 
READ (9, 112)  IBLA,DELPA, SBA, EBA 

111  FORMAT  (2I5,3E10. 3) 

112  FORMAT  (15, 3E10. 3) 

115  FORMAT  (4(15, E10.  3) ) 

NPLT=8 

DO  10  ID=1,NPLT 

READ(g, 120)  (ITL(I) , 1=1,8) 

120  FORMAT (aA10) 

121  FORMAT (4A10) 

READO,  121)  (LABX  (I) ,  1  =  1,  4) 

READ(g, 121)  (LABY(I), 1=1,4) 

READO,  130)  NEPTS,  XPMN,  XPMX,YPMN,  YPMX 
130  FORMAT (15, 4E15.  8) 

READO,  140)  (XARY  (K) ,  K=l,  NEPTS) 

READ (9, 140)  (YARY(L).L=1,  NEPTS) 

140  FORMAT (10E 15. 8) 

IK=1 

IF  (ID. LT.6)  GO  TO  10 

IF  (ID.  LE.2.  OR.  ID.  EQ.  6.0R.  ID.  EQ.  8)  GO  TO  151 
C  NEPTS=150 

C  XPMX=2990. 

151  CONTINUE 

CALL  PLOTTER (XARY,  YARY,  NEPTS,  XPMN, XPMX,  YPMN,  YPMX,  IK,  LABX,  LABY) 
CALL  ENDPL(-l) 

10  CONTINUE 
CALL  GDONE 
END 
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SUBROUTINE  PLOTTER < XflRY, YftRY, NP, XMN, XMX, YMN,  YMX, KIND, LPBX, LfiBY) 
C 

COMMON  /PLOTV/  ITL (8) , ISTL (8) , IDB 
COMMON  /PIT)/  NPTS,  NSKIP,  TFAC,  XFflC,  «FflC, 

*ISPBX,  TSPBX,  ISPEX,  TSPEi^  ISPBA,  TSPBA,  ISPEA,  TSPEA, 

*IBLX, DELPX, SBX, EBX,  IBLA,  DELPA,  SBA,  EBA 
DIMENSION  XARY (NP) , YARY (NP) ,  LABX (A) ,  LABY ( A) 

WRITE (6, £300)  NP,  XMN,  XMX,  YMN, YMX, KIND 
£300  FORMAT (5X,*  ENTERED  PLOTTER 

*  *NP,  XMN,  XMX, YMN, YMX, KIND  =  *,  15, A ( 1 X, E10. A) , 15) 

CALL  HEIGHT (0.1) 

IF(KIND. LT. 0)  GO  TO  £00 
C 

IF(KIND. EQ. £)  GO  TO  100 

*****  IF  KIND. EQ. 1  THEN  PLOT  IS  LINEAR-LINEAR  ***** 

50  LINET  =  0 
LINES  =  0 

CALL  SCLl (XMN,  XMX,  XORG,  XSTP,  XEND) 

CALL  SCLl (YMN,  YMX, YORG,  YSTP,  YEND) 

WRITE (6, £303)  XORG, XSTP, XEND, YORG, YSTP, YEND 
S303  FORMAT (£X, ^LINEAR  PLOT  *, 6 <£X, E15. 6) ) 

CALL  RLINER ( XORG, XSTP, XEND, YORG, YSTP, YEND, LABX, LABY) 

CALL  DRAWC (XARY, YARY, NP, LINET, LINES) 

GO  TO  A00 

*****  IF  KIND. EQ. £  THEN  PLOT  IS  LOG-LOG  ***** 

100  LINET  =  0 
LINES  =  0 

CALL  l>CL£  ( XMN,  XMX,  XORG,  XCYC,  KIND) 

IF(KIND. EQ. 1)  GO  TO  50 
CALL  SCL£ ( YMN, YMX, YORG,  YCYC,  KIND) 

IF(KIND. EQ. 1)  GO  TO  50 
WRITE (6, £305)  XORG, XCYC, YORG, YCYC 
£305  FORMAT (5X, *L06-L0G  PLOT  *,  A (£X,  E15.  6) ) 

CALL  LOGLLL ( XORG, XCYC, YORG, YCYC, LABX, LABY) 

CALL  DRAWC (XARY, YARY, NP, LINET, LINES) 

GO  TO  A00 

*****  IF  KIND.LT. 0  THEN  PLOT  AN  OVERLAY  ***** 

£00  LINET  =  LINET+1 
WRITE (6, £307) 

£307  FORMAT (5X,*  OVERLAY  PLOT  *) 

CALL  BLOFF(IDB) 

CALL  MESSAGdSTL,  80,  0.  0,  6.  £5) 
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CPILL  MESSfiG<4HDPlTft,  4,4.5,  5.  8) 

COLL  STRTPT<5.S,  5.  8) 

COLL  C0NNPT<5.  8,5.  8) 

COLL  MESS0G(4HFIT  ,4,4.5,  5.  6) 

COLL  DOSH 

COLL  STRTPT(5.  £,5.6) 

COLL  C0NNPT<5.  8,5.6) 

COLL  RESET (4HD0SH) 

COLL  DROWC  (XORY,  YORY,  r4P,  LINET,  LINES) 

400  CONTINUE 

COLL  MESSOG(30HNO.  OF  PTS  REOD  FROM  TOPE 
COLL  INTNO(NPTS, 8. a,  5.  5) 

COLL  MESSOG(30HSKIP  FOCTOR 
COLL  INTNQ(NSKIP,  a.  a,5.£5) 

COLL  MESSOG<30HTIME  CONVERSION  FOCTOR 
COLL  REOLNOCTFOC,  a,  8.8,5.  ) 

COLL  MESSOG <30H INPUT  DOTO  CONV.  FOCTOR 
COLL  REOLNO  <  XFOC,  2,  8.  8,  4.  75) 

COLL  MESSOG <30HOUTPUT  DOTO  CONV.  FOCTOR 
COLL  REOLNO<OFOC,  2,  a.  8,  4.5) 

TNP=4. 5 

IF  <ISPBX. NE. 1)  GO  TO  11 
TNP=TNP-. 25 

COLL  MESSOG <30HEEG.  INPUT  SPLINE  BEFORE 
COLL  RE0LN0<TSPBX,4,8.a,TNP) 

11  IF  <ISPEX.NE. 1)  GO  TO  12 
TNP=TNP-. £5 

COLL  MESSOG <30HEND  INPUT  SPLINE  OFTER 
COLL  RE0LN0(TSPEX,4,a.  S,TNP) 

12  IF  (ISPBO. NE. 1)  GO  TO  13 
TNP=TNP-. 25 

COLL  MESSOG <30HBEG.  OUTPUT  SPLINE  BEFORE 
COLL  REOLNO (TSPBO,  4,  8.8, TNP) 

13  IF  <ISPEO. NE. 1)  GO  TO  14 
TNP=TNP-. £5 

COLL  MESSOG <30HEND  OUTPUT  SPLINE  OFTER 
COLL  REOLNO<TSPEO,  4,  8.  a,  TNP) 

14  IF  (IBLX.NE. 1)  GO  TO  15 
TNP=TNP-. 25 


,  30,  6 
,  30,  6 
,  30,  6 

,  h^0,  6 

,  30,  6, 


,  30,6 


,  30,  6. 


,  30,  6. 


,  30,  6. 


COLL  MESSOG <30HINPUT  DOTO  BOSELINE  CORRECTION, 30, 6. 


TNP=TNP-. £5 

COLL  MESSOG <30H 

COLL  REOLNO(DELPX,£,  a.  8,TNP) 

TNP=TNP-. £5 

COLL  MESSOG (30H 

COLL  RE0LL'0<SBX,4,  a.  8,  TNP) 

TNP=TNP-. 25 

COLL  MESSOG (30H 

COLL  RE0LN0(EBX,4,a.  8,  TNP) 


DELPX 

=  ,30,6 

SBX 

=  ,30,6 

EBX 

=  ,  30,  6 

1,5.  5) 

1.5.  £5) 

1.5. ) 
1,4.  75) 
1,4.  5) 

1, TNP) 

1, TNP) 

1,TNP) 

1, TNP) 

1, TNP) 
1, TNP) 

1, TNP) 

1, TNP) 
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15  IF  (IBLfl. NE, 1)  GO  TO  900 
TNP=TNP-. £5 

COLL  MESSflGOlHOUTPUT  DflTFl  BASELINE  CORRECTION,  31, 6.  1,  TNP) 
TNP=TNP-. £5 

CALL  MESSAG(30H  DEt.PA  =  ,30, 6.1,  TNP) 

CALL  REALNO(DELPt  ,£,  8.  8,TNP) 

TNP=TNP-.  £5 

CALL  MESSAG(30H  SBA  =  ,30, 6.1, TNP) 

CALL  REALNO(SBA, 4,8.  8,  TNP) 

TNP=TNP-.  £5 

CALL  MESSAG(30H  EBA  =  ,30, 6.1, TNP) 

CALL  REALNO (EBA, 4, 8.  8,  TNP) 

900  CONTINUE 
RETURN 
END 

SUBROUTINE  SCLKXMN,  XMX,  AORB,  ASTP,  AMAX ) 

DIMENSION  S(7) 

C 

C  *****  FIND  LINEAR  SCALES  ***** 

C 

WRITE (6,  £300)  XMN, XMX 

£300  FORMAT <5X,*SUBR0UTINE  SCLl  XMN, XMX  =  *,  £ (E15. 6, £X) ) 

SMIN  =  0.00006 
S(l)  =  0. 0001£ 

S<£)  «  0.00018 
S<3)  =  0.  000£4 
S(4)  =  0.00030 
S(5)  =  0.00036 
S(6)  =  0.00060 
S(7)  =  0.  001£0 
C 

DIF  =  XMX  -  XMN 
IF(DIF.LT.S(1) )  GO  TO  90 
5  CONTINUE 
DO  10  1  =  1,7 
lU  =  I 

10  IF(DIF.LT.S(I)  )  GO  TO  30 
DO  £0  J=l,7 
£0  S(J)  =  S(J)*10.0 

IF(S<1).GT. 1.0E15)  STOPlll 
GO  TO  5 
C 

30  DMAX  =  S(IU) 

DSTP  =  DMAX/6. 0 
C 

C  DETERMINE  OFFSET 

C 

IF(XMN. LT. 0. 0)  GO  TO  60 
DORG  =  0.  0 

IF (XMN. LT. DSTP)  GO  TO  99 
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OFFSET  =  DSTP 
35  OFFSET  =  OFFSET+DSTP 

IF (XMN.6T. OFFSET)  GO  TO  35 
DORG  =  OFFSET-DSTP 
DMftX  =  DMftX+DORG 
GO  TO  99 

60  OFFSET  =  0.0 
65  OFFSET  =  OFFSET-DSTP 

IF(XMN.LT. OFFSET)  GO  TO  65 
DORG  =  OFFSET 
DMftX  =  DMftX+DORG 
IF (XMX. LT. DMftX)  GO  TO  99 
IFdU.LT.?)  DMftX  =  S(IU+1) 

IF<IU. EQ. 7)  DMftX  =  S(2)*10. 0 
DSTP  =  DMftX /6.0 
GO  TO  60 

DIFFERENCE  IS  ZERO 

90  CONTINUE 

DORG  =  XMN-SMIN 
DMftX  =  XMN+SMIN 
DSTP  =  SMIN/3. 0 

99  ftORG  =  DORG 
ftSTP  =  DSTP 
ftMftX  =  DMftX 

WRITE (6, 2303)  DORG, DSTP, DMftX 
2303  F0RMftT(5X,*  LEftVING  SCLl  *,  3 (E15.  6,  EX) ) 

RETURN 
END 

SUBROUTINE  SCL2(XMN,  XMX,  ftORG, ftCYC, KIND) 

SCftLE  FOR  LOG-LOG  PLOTS 

WRITE (6,  2300)  XMN,  XMX 
2300  F0RMftT(5X,*ENTER  SCL2  *, 2 (E15. 6, EX ) ) 

IFCXMN. LT.  1.0E-a)  GO  TO  80 
IF(XMX.LT.  1.0E-8)  GO  TO  81 
C 

SMN  =  HLOG10(XMN) 

SMX  =  flLOG10<XMX) 

MN  =  IFIX<SMN) 

IF<SMN. LT.0. 0)  MN=MN-1 
MX  =  IFIX(SMX) 
ftORG  =  10. **MN 
DIF  =  (MX-MN)+1 

IF<MN. LT. 0  .ftND.  MX. LT. 0)  DIF  =  MX-MN 
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ttCYC  =  ftBS<6.0/DIF) 
GO  TO  90 


80  WRITE <6, 1000)  XMN 


1000 

FORMAT (5X,*XMN  =  *, 
BO  TO  8£ 

E15.  6,  * 

A 

LINEAR 

PLOT 

WILL 

BE 

MADE. *) 

81 

WRITE (6, 1001)  XMX 

1001 

FORMAT (5X, *XMX  =  *, 

E15.  6,  * 

A 

LINEAR 

PLOT 

WILL 

BE 

MADE. *) 

8£  KIND  =  1 
C 


90  CONTINUE 

WRITE  <6,  £303)  MN,  MX,  DIF, ftORG, flCYC 

£303  FORMAT <5X,*LEAVING  SCL£  MN,  MX,  DIF,  AORG,  ACYC*,  £15,  3 < IX,  E15. 6) ) 
RETURN 
END 

SUBROUTINE  DRAWC (X,  Y,  NP,  LINET,  LINES) 

DIMENSION  X<NP),Y<NP) 

C 

WRITE (6, £300)  NP,  LINET,  LINES 
£300  FORMAT (5X,  CENTER  DRAWC  NP, LINET, LINES  =  *,315) 

IF (LINET. LE. 0)  GO  TO  10 
IF(LINET.EQ- 1)  CALL  DASh 
IF (LINET. EQ. £)  CALu  CHNDOT 
IF (LINET. EQ. 3)  CALL  CHNDSH 
IF (LINET. EQ. 4)  CALL  DOT 
C 

10  CALL  CURVE(X,Y,NP, LINES) 

C 

IF (LINET. LE. 0)  GO  TO  99 
CALL  RESET (3HALL) 

CALL  HEIGHT (0.1) 

C 

99  CONTINUE 
RETURN 
END 

SUBROUTINE  RLINER (XORG, XSTP, XEND,  YORG, YSTP,  YEND, LABX, LABY) 

C  COMMON  /COUNT/  ICOUNT, lOPT, LFILT 

COMMON  /PLOTV/  ITL (8) ,  ISTL (8) , IDB 
DIMENSION  LABX (3), LABY (3) 

C 

WRITE (6, £300) 

£300  FORMAT  (5X,  *ENTERED  RLINER . *) 

CALL  PAGE (10. 5, 8. 5) 

CALL  PHYSOR(1.0,  1.0) 

CALL  XNAME(LABX, 30) 

CALL  YNAME(LABY,30) 

CALL  AREA£D(6.  0,6.  0) 

IFdOPT.  EQ.  1)  CALL  BLREC  (4.  4,  5.  5,  1 . 6,  0.  5,  1 . 0) 

I F ( I OPT. EQ. 1 )  CALL  BLKEY (IDB) 

CALL  MESSAG(ITL,80, 0. 0, 6. 5) 
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CPILL  GRAF  ( XORG,  XSTP,  XEND,  YORG,  YSTP,  YEND) 

CALL  DOT 
CALL  GRID <1,  1) 

CALL  RESET (3HD0r> 

C 

RETURN 

END 

SUBROUTINE  LOGLLL(XOR,  XCY, YOR, YCY,  LABX,  LABY) 

C  COMMON  /COUNT/  ICOl’NT,  lOPT,  LFILT 

COMMON  /PLOTV/  ITL <8) , ISTL (fl) , IDB 
DIMENSION  LABX (3) , LABY (3> 

C 

WRITE <8, 2300) 

2300  FORMAT  <5X,  ^ENTERED  LOGLLL . 

CALL  PAGE ( 10.  5,  8.  D) 

CALL  PHYSORd.  0,  1.0) 

CALL  XNAME<LABX,30) 

CALL  YNAME<LABY,30> 

CALL  AREA2D(6.  0,6.0) 

C  IFdOPT.EQ.  1)  CALL  BLPEC<4.  4,  5.  5,  1. 6,  0.  5,  1^  0) 

C  IFdOPT.EQ.  1)  CALL  BLKEYdDB) 

call.  MESSAG<ITL,  80,0.0,  6.  5) 

CALL  LOGLOG<XOR,  XCY,  YOR,  YCY) 

CALL  DOT 
CALL  GRID  < 1,1) 

CALL  RESET (SnDGT) 

C 

RETURN 

END 
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Input  x(t) 


I - - - 

Linear  Transfer  Function 


^  Output 


Figure  1.  Linear  cause-effect  relationship  between  two  data  records. 


Pressure 


Time  (ms)  Frequency  (kHz) 

(a)  Record  #6  time  history  (b)  Record  #6  FFT 


Time  (ms) 


(c)  Record  #6  after  direct  Inverse  FFT 


Figure  3.  Illustration  of  an  original  time  history,  an  FFT,  and  an  Inverse  FFT 


Impulse  (psi-s) 


Figure  5.  Record  #2,  REST  pressure  on  the  structure. 
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Impulse  (psi-s) 


'x>  ^  ff  ^.^'ik.. 


Figure  6. 


Record  #4, 


HEST  pressure  on  the  soil. 
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Figure  g.  Record  #7,  soil  pressure  at  5.21'  depth. 
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Figure  11.  Record  §9,  vertical  structure  acceleration  at  3.28'  depth. 
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Figure  13.  Record  #11,  vertical  soil  accele'^ation  at  12.21'  depth. 
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'0  5  10  15  20  25  30  35  40 

Time  (ms) 

Figure  19.  Record  #17,  vertical  structural  strain  at  18 


90 


RfiNGE  (KM)  -  0.04554 

POS,  PHRSE  I  SEC)  -  0.35691 

TOR  (SEC)  -  0.00134 

LOW  PASS  no  (HZ)  -  200. 
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Figure  23b.  Time  history  comparison  of  Speicher-Brode  waveform  to  REST  record  #2. 
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RfiNGE  (KM)  -  0,03772 
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I  tffU 


Figure  24a.  FFT  fit  comparison  for  test  ciata  record  #4. 


PRESSURE  HISTORY 


waveform  to  HEST  record  #4 


rOURFIT  SPEICHER  ERODE 


Figure  24c.  Impulse  time  history  comparison  for  record  #4 
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Figure  25.  FRF's  for  vertical  soil  s 
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Figure  29.  FRF's  for  structural  hoop  strains  (1.3',  5.3',  and  12.2 
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Figure  30,  Comparison  of  orio 
record  #6. 
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Figure  33.  Program  FREQRES  flow  chart  (Continued) 
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Spline  End 
Portion  of  th 
Output  Data 
Record 
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Figure  33.  Program  FI^EQRES  flow  chart  (Continued) 
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Save  inverse 
fFFT  of  Output 
on  Plot  File 
(Tape48) 


FFT  the 

ppeicher-Brode 
Fit  to  Input 
Data 


Multiply  FRF 
times  Speicher-j 
prode  FFT  to 
get  Modified 
flutpiit.-EE'T _ 


>ave  Modified 
(Output  FFT  on 
Plot  Filej 
(Tape48) 


Inverse  FFT 
[lodified  Output 
FFT  to  get  nev 
putput  time 

h-is.to 


Save  new  Out- / 
(put  time  his¬ 
tory  on  Plot 
/File  (Tape48) 


Figure  33.  Program  FREQRES  flow  chart  (Concluded). 


Table  1.  Descriptions  of  Example  Test  Data 


Record  Number 

2 

4 


Description 

HEST  pressure  on  the  test  article 
REST  pressure  on  soil 

stress  at  0,5'  depth 
stress  at  0.5’  depth 
stress  at  5.21'  depth 
structure  acceleration  at  0.83'  depth 
structure  acceleration  at  3.28'  depth 
soil  acceleration  at  5.21'  depth 
soil  acceleration  at  12.21'  depth 


soil 

soil 


Vertical 
Vertical 
Vertical  soil 
Vertical 
Vertical 
Vertical 
Vertical 
Vertical 
Vertical 
Vertical 
Vertical 
Vertical 
Vertical 


structure  strain  at  1.29' 
structure  strain  at  1.29' 
structure  strain  at  3.33' 
structure  strain  at  5.21' 
structure  strain  at  12.21' 

.  structure  strain  at  18.71'  depth 
Structure  hoop  strain  at  1.29'  depth 
Structure  hoop  strain  at  5.29'  depth 
Structure  hoop  strain  at  12.21'  depth 


depth 

depth 

depth 

depth 

depth 


Table  2.  FREQRES  input  file  variable  list  and  descriptions 


wARD 

COLUMN 

FORMAT 

VARIABLE 

DESCRIPTION 

1 

1-5 

IS 

NEPTS 

NO.  OF  POINTS  TO  BE  READ  FROM  TAPE 

6-10 

15 

NSKIO 

SKIP  INTERVAL  (DEFAULT'1.  ALL  POINTS 
FROM  TAPE  ARE  SAVED) 

11-20 

E10.3 

TFAC 

TIME  CONVERSION  FACTOR  (DEF AULT- 1 . 0) 

21-30 

E10.3 

XFAC 

INPUT  DATA  CONV.  FACTOR  ( DE F AULT- 1 . 0) 

31-40 

E10.3 

AFAC 

OUTPUT  DATA  CONV.  FACTOR  ( DEFAULT- 1 . 0 ) 

2 

1-5 

15 

ISPBX 

0:  NO  SPLINE  PERFORMED  ON  BEGINNING 

Qc  input  data 

1:  BEGINNING  OF  INPUT  DATA  WILL  BE 
SPLINED  (DEFAULT-0) 

6-15 

E1C.3 

TSPBX 

IF  ISPBX- 1.  A  SPLINE  IS  PERFORMED  AT 
THIS  TIME  BACK  TO  TIME  ZERO  (DEFAULT 

IS  0.0  WHICH  MEANS  NO  SPLINE  IS  DONE) 

16-20 

15 

ISPEX 

0  NO  SPLINE  PERFORMED  ON  END  OF 

INPUT  DATA 

1-  END  DF  INPUT  DATA  WILL  BE  SPLINED 
(DEFiuLT-O) 

21-30 

E10.3 

TSPEX 

IF  ISPEX- 1.  TIME  AT  WHIC:i  SPLINE 

BEGINS  FOR  INPUT  DATA  (DEFAULT  IS 

85%  OF  TTOT) 

31-35 

15 

ISPBA 

0:  NO  SPLINE  PERFORMED  ON  BEGINNING 

OF  OUTPUT  DATA 

1:  BEGINNING  OF  OUTPUT  DATA  WILL  BE 
SPLINED  (DEFAULT-0) 

36-45 

E10.3 

TSPBA 

IF  lSPBA-1,  A  SPLINE  IS  PERFORMED  AT 
THIS  TIME  BACK  TO  TIME  ZERO  (DEFAULT 

IS  0.0  WHICH  MEANS  NO  SPLINE  IS  DONE) 

46-50 

15 

ISPEA 

O:  NO  SPLINE  PERFORMED  ON  END  OF 

OUTPUT  DATA 

Is  END  OF  OUTPUT  DATA  WILL  BE 

SPLINED  (DEFAULT-0) 

51-60 

E10.3 

TSPEA 

IF  ISPEA-1.  TIME  AT  WHICH  SPLINE 

BEGINS  FOR  OUTPUT  DATA  (DEFAULT 

IS  85%  OF  TTOT) 

3 

1-5 

15 

IBLX 

INPUT  DATA  BASELINE  CORRECTION  TRIGGER 
0:  NO  BASELINE  CORRECTION 

1:  BASELINE  CORRECTION  WITH  THE 

FOLLOWING  PARAMETERS  (DEFAULT-0) 

6-15 

E10.3 

DELPX 

CORRECTION  ADDED  TO  INPUT  DATA  VALUES 
AFTER  TIME  SBX.  IF  EBX  AND  SBX  ARE 
EQUAL  THEN  THE  FULL  VALUE  OF  DELPX  IS 
AODED  AT  ALL  TIMES  AFTER  SBX.  IF  EBX 

IS  GREATER  THAN  SBX  THEN  THE  PLOT  IS 
ROTATED  ABOUT  THE  POINT  DEFINED  AT  SBX 
BY  THE  AMOUNT  DEI PX  AT  TIME  EBX. 

16-25 

E10.3 

SBX 

start  time  for  BASELINE  CORRECTION 
(PLOT  ROTATION  POINT  IF  EBX  >  S3X) 

26-35 

E10.3 

EBX 

END  TIME  FOR  BASELINE  CORRECTION 

36-40 

15 

INTX 

INPUT  DATA  INTEGRATION  TRIGGER 

0:  NO  INTEGRATION 
1:  INTEGRATE  INPUT  DATA 


Table  2.  FREQRES  input  file  variable  list  and  descriptions  (Concluded). 


4 

1-5 

15 

I8LA 

OUTPUT  data  baseline  CORRECTION  TRIGGER 
0:  NO  BASELINE  CORRECTION 

1:  BASELINE  CORRECTION  WITH  THE 

rCLLOWlNG  I'ARAMETERS  {DEFAULTED) 

6-15 

E10,3 

DELPA 

CORRECTION  ADDED  TO  OUTPUT  DATA  VALUES 
AFTER  TIME  SBA.  IF  EBA  AND  SBA  ARE 

EQUAL  THEN  THE  FULL  VALUE  OF  DELPA  IS 
ADDED  AT  ALL  TIMES  AFTER  SBA.  IF  EBA 

IS  GREATER  THAN  SBA  THEN  THE  PLOT  IS 
ROTATED  ABOUT  THE  POINT  DEFINED  AT  SBA 
BY  TV!E  amount  DELPA  AT  TIME  EBA. 

16-25 

E10.3 

SBA 

START  TIME  FOR  BASELINE  CORRECTION 
(PLOT  ROTATION  POINT  IF  EBA  >  SBA) 

26-35 

E10.3 

EBA 

END  TIME  FOR  BASELINE  CORRECTION 

36-40 

15 

INTA 

OUTPUT  DATA  INTEGRATION  TRIGGER 

0:  NO  INTEGRATION 

1:  INTEGRATE  OUTPUT  DATA 

NOTE: 

ALL  OF 

THE  FOLLOWING 

LABELS  SHOULD  BE  CENTERED  WITHIN  THE 

FIRST 

30  COLUHNS  OF  EACH  LINE  OF  THE  INPUT  FILE. 

5 

1-40 

4A10 

ITX 

INPUT  DATA  X-AXIS  cABEL:  EXAMPLE: 

TIME  (SEC) 

6 

1-40 

4A10 

ITY 

INPUT  DATA  Y-AXlS  LABEL:  EX: 

PRESSURE  (PSD 

7 

1-40 

4A10 

ITX 

OUTPUT  DATA  X-AXIS  LABEL:  EX: 

TIME  (SEC) 

e 

1-40 

4A10 

ITY 

OUTPUT  DATA  Y-AXIS  LABEL:,  EX: 

STRAIN  (IN/IN) 

9 

1-40 

4A  10 

ITX 

INPUT  DATA  FOURIER  AMP.  SPEC. 

X-AXIS  LABEL:  EX: 

FREQUENCY  (HZ) 

10 

1-40 

4A10 

ITY 

INPUT  DATA  FOURIER  AMP.  SPEC. 

Y-AXIS  LABEL:  ex- 

amplitude  (PSI-SEC) 

1 1 

1-40 

4A10 

ITX 

OUTPUT  DATA  FOURIER  AMP.  SPEC. 

X-AXIS  LABEL:  EX - 

FREQUENCY  (HZ) 

12 

1-40 

4A10 

ITY 

OUTPUT  DATA  FOURIER  AMP.  SPEC. 

Y-AXIS  LABEL:  EX: 

amplitude  (SEC) 

13 

1-40 

4A10 

ITX 

FREQUENCY  RESPONSE  FUNCTION 

X-AXIS  LABEL:  EX: 

FREQUENCY  (HZ) 

14 

1-40 

4A10 

ITY 

FREQUENCY  RESPONSE  FUNCTION 

Y-AXIS  LABEL:  EX: 

FFTX/FFTA 

15 

1-40 

4A10 

ITX 

ERODE  OUTPUT  RESPONSE 

X-AXIS  LABEL:  EX. 

TIME  (SEC) 

16 

1-40 

4A10 

ITY 

BRODE  OUTPUT  RESPONSE 

Y-AXIS  LABEL:  ex¬ 
strain  (IN/IN) 


Table  3 


Sample  output  listing  from  a  FREQRES  calculation 


FREQRES  OUTPUT  LISTING 

THE  NUMBER  OF  POINTS  READ  FROM  THE  DATA  RECORO  TAPES  IS  9370  WITH  A  SKIP  OF  27  CONSIDERED  FOR  ANALYSIS 

TIME  CONVERSION  FACTOR  =  . 100E+01 

INPUT  DATA  CONV.  FACTOR  =  . lOOE+Oi 
OUTPUT  DATA  CONV.  FACTOR  =  . 100E+01 

BEGINNING  PORTION  OF  INPUT  DATA  SPLINED  FROM  TIME  EQUAL  0.0  TO  TIME  EQUAL  .260E-02 
FINAL  15%  OF  INPUT  DATA  SPLINED  TO  ZERO 

BEGINNING  PORTION  OF  OUTPUT  DATA  SPLINED  FROM  TIME  EQUAL  0.0  TO  TIME  EQUAL  .300E-02 
FINAL  15%  OF  OUTPUT  DATA  SPLINED  TO  ZERO 


********************************************** 


VERY  IMPORTANT  NOTICE:- 

NBPTS  FOR  THE  SPEICHER-BRODE  FIT  FROM  FOURFIT  =  392 

MUST  BE  EQUAL  TO  NPT  FROM  THIS  PROGRAM  =  369 

IF  NBPTS  AND  NPT  ARE  NOT  EQUAL.  THIS  PROGRAM  WILL 
TRUNCATE  ONE  OF  THEM  TO  MAKE  THEM  EQUAL. 

OTBP  FOR  THE  SPEICHER-BRODE  FIT  FROM  FOURFIT  =  . 13372 120E -03 

MUST  BE  VERY  CLOSE  TO  DTD  FROM  THIS  PROGRAM  =  . 13500000E-03 

IF  THEY  ARE  NOT  CLOSE,  THIS  PROGRAM  WILL  STOP 


OUTPUT  RECORD  TOTAL  IMPULSE 
FIRST  POINT  OF  OUTPUT  FFT 
OUTPUT  INVERSE  FFT  OFFSET 


. 64978 1E+02 
.648020E+02 
.  129430E+04 
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Base  Line  Correction 


Table  5.  Summary  of  the  total  impulse,  first  FFT  value,  and 
inverse  FFT  offset  for  the  test  data. 


Record  Total  Impulse  First  FFT  Value  Inverse  FFT  Offset 

Number 


4 

24.93  psi-sec 

24.87  psi-sec 

516  psi 

5 

64.98  psi-sec 

64.80  psi-sec 

1294  psi 

6 

47.34  psi-sec 

47.22  psi-sec 

947  psi 

7 

33.74  psi-sec 

33.65  psi-sec 

676  psi 

8 

.5192 

ft 

.5174 

ft 

10.4  ft/s 

9 

.7249 

ft 

.7223 

ft 

14.5  ft/s 

10 

.7842 

ft 

.7821 

ft 

15.7  ft/s 

11 

.4031 

ft 

.4020 

ft 

8.07  ft/s 

12 

.1170 

mi  cro 

sec 

.1166 

micro 

sec 

-.39  micro 

strai n 

13 

.6474 

mi  cro 

sec 

.6451 

micro 

sec 

10.2  micro 

strain 

14 

-14.46 

mi  cro 

sec 

14.40 

micro 

sec 

-288  micro 

strai n 

15 

8.581 

micro 

sec 

8.551 

micro 

sec 

173  micro 

strain 

16 

48.47 

micro 

sec 

48.30 

micro 

sec 

968  micro 

strain 

17 

5.014 

micro 

sec 

4.997 

micro 

sec 

98.4  micro 

strain 

18 

24.80 

micro 

sec 

24.71 

micro 

sec 

496  micro 

strain 

19 

10.27 

mi  cro 

sec 

10.23 

micro 

sec 

205  micro 

strain 

20 

-4.586 

micro 

sec 

4.570 

micro 

sec 

-92.1  micro 

strain 

Table  6. 


FREQRES  input  files. 


Description 


Input  file  for  input  variable 
specification  as  described  in 
Section  3.1 

Spei cher-Brode  "best- fit" 
waveform  obtained  from  TAPE49 
of  FOURFIT  calculation 

Input  data  record  digitized 
time  history 

Output  data  record  digitized 
time  history 


FREQRES  output  files. 


Description 
Printed  output 
Plot  file 


DfSTRIBUTION  LIST 


DEPARTMENT  OF  DEFENSE 

DEFENSE  INTELLIGENCF  AGENCY 
ATTN  DB-4C 
ATTN  RTS-2A  TECH  LIB 
ATTN.  RTS-2B 

DEFENSE  NUCLEAR  AGENCY 
2  CYS  ATTN  SPSS 
4  CYS  ATTN.  STTI-CA 

DEFENSE  TECHNICAL  INFORMATION  CENTER 
12  CYS  ATTN  DO 

DEPARTMENT  OF  DEFENSE  tXHLO  SAFETY  BOARD 
ATTN  CHAIRMA.'I 

FIELD  COMMAND,  DNA.  UET  2 
LAWRENCE  LI'.'ERMORE  NATIONAL  LAB 

Ar:j  rc-1 

FirXD  COMMAND,  DNA 
ATTN.  FCPR 
ATTN  FCT 

ATTN:  FCTTWSUMMA 
ATTN  FCTXE 

JOINT  STRAT  TGT  PLANNING  STAFF 
ATTN  JLK  (ATTN.  DNA  REP) 

UNDER  SECY  OF  DEF  FOR  RSCH  &  ENGRG 
ATTN-  STRAT  &  SPACE  SYS(OS) 

DEPARTMENT  OF  THE  ARMY 

BMD  ADVANCED  TECHNOLOGY  CENTER 
ATTN.  ATC-T 
AITN  ICRD.ABH-X 

BMD  SYSTEMS  COMMAND 

ATTN:  BMOSC-LEE  R  BRADSHAW 
ATTN  BMDSC-LEH  R  C  WEBB 

DEP  CH  OF  STAFF  FOR  OPS  &  PLANS 

ATTN  DAMO-NCNUCCHEMDIR 

HARRY  DIAMOND  UBORATORIES 

ATTN.  DELHD-TA-L  81100,  TECH  LIB 
ATTN  SLCHD-NW-P 

U  S  ARMY  ARMAMENT  MATERcAL  READINESS  COMMAND 
ATTN  MA  LIBRARY 

U  S  ARMY  BALLISTIC  RESEARCH  LAB 
2  CYS  ATTN  DRDAR-BLA-S  TECH  LIB 
ATTN.  DRDAR-BLT  J  KEEFER 


U  S  ARMY  ENGINEER  DIV  HUNTSVILLE 
ATTN:  HNDED-SR 
ATTN:  HNDEN-FO 

U  S  ARMY  ENGR  WATERWAYS  EXPER  STATION 
ATTN:  J  STRANGE 
ATTN:  JZEUSKO 
ATTN.  LIBRARY 
ATTN  WESSDJ  JACKSON 
ATTN  WESSE 

U  S  ARMY  FOREIGN  SCIENCE  &  TECH  CTR 
ATTN.  DRXST-SD 

U  S  ARMY  NUCLEAR  &  CHEMICAL  AGENCY 
ATTN  LIBRARY 

US  ARMY  WHITE  SANDS  MISSILE  RANGE 
ATTN  STEWSTE-N  K  CUMMINGS 

USA  MISSILE  COMMAND 

ATTN-  DOCUMENTS  SECT  ION 
ATTN:  DRSMI-RH 

DEPARTMENT  OF  THE  NAVY 

DAVID  TAYLOR  NAVAL  SHIP  R  &  D  CTR 

ATTN  TECH  INFO  CTR  CODE  522  1 

NAVAL  CIVIL  ENGINEERING  LABORATORY 
ATTN:  CODEL-51 

NAVAL  FACILITIES  ENGINEERING  COMMAND 
ATTN  CODE04B 

NAVAL  POSTGRADUATE  SCHOOL 

ATTN:  CODE  1424  LIBRARY 

N.AVAL  RESEARCH  LABORATORY 

ATTN:  CODE  2627  TECH  LIB 

NAVAL  SURFACE  WEAPONS  CENTER 
ATTN:  C0DEF31 

NAVAL  SURFACE  WEAPONS  CENTER 

ATTN  tech  LIB  &  INFO  SVCS  BR 

NAVAL  WF.APONS  CENTER 

ATfN.  CODE  266  C  AUSTIN 
ATTN.  CODE  343  FKA6A2,  TECH  SVCS 

NAVAL  WE.AI’ONS  EVALUATION  FACILITY 
ATTN.  CLASSIFIED  LIBRARY 

OFFICE  OF  NAVAL  RESEARCH 
ATTN.  CODE  474 


U  S  ARMY  CORPS  OF  ENGINEERS  DEPARTMENT  OF  THE  AIR  FORCE 

ATTN.  DAcN-ROL 

AFRCE-BMS/DEE 

U  S  ARMY  ENGINEER  DIV  OHIO  RIVER  ATTN  DEB 

ATTN  ORDAS  L  TECH  LIB 


DEPARTMENT  OF  THE  AIR  FORCE  (CONTINUED) 

AIR  FORCE  GEOPHYSICS  UBORATORY 
ATTN;  LWH  MOSSING 

AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 
ATTN  LIBRARY 

AIR  FORCE  SYSTEMS  COMMAND 
ATTN.  DLW 

AIR  FORCE  WEAPONS  LABORATORY,  AFSC 
ATTN-  NTE  M  PLAMONDON 
ATTN:  NTED  R  MATALUCCI 
ATTN.  NTES-CRHENNY 
ATTN.  SUL 

BALLISTIC  MISSILE  OFFICE/DAA 
ATTN.  ENSN 

ATTN.  MGEN  A  SCHENKER 
ATTN  PP 

DEPUTY  CHIEF  OF  STAFF 
ATTN:  AF/RDQI 

FOREIGN  TECHNOLOGY  DIVISION,  AFSC 
ATTN  Nils  LIBRARY 

STRATEGIC  AIR  COMMAND 
ATTN.  NRI/STINFO 

STRATEGIC  AIR  COMMAND 
ATTN  XPFS 

DEPARTMENT  OF  ENERGY 

DEPARTMENT  OF  ENERGY 

ALBUQUERQUE  OPERATIONS  OFFICE 
ATTN  CTID 
ATTN  R  JONES 

DEPARTMENT  OF  ENERGY 

OFC  OF  MIL  APPLICATION,  GTN 
ATTN  OMA/RD&T 

DEPARTMENT  OF  ENERGY 

NEVADA  OPERATIONS  OFFICE 

ATTN  DOC  CON  FOR  TECH  LIB 

UNIVERSITY  OF  CALIFORNIA 

LAWRENCE  LIVERMORE  NATIONAL  LAB 
ATTN-  L-203RSCHOCK 
ATTN  L-203  T  BUTKOVICH 
ATTN  L-658  TECH  INFO  DEPT  LIB 
ATTN.  P COYLE 

LOS  ALAMOS  NATIONAL  LABORATORY 
ATTN.  MSP364RPTSUB 
ATTN  MS530  G  SPILLMAN 
ATTN  R  WHITAKER 

OAK  RIDGE  NATIONAL  LABORATORY 

ATTN  CENTRAL  RSCH  LIBRARY 
ATTN  CIVIL  DEF  RES  PROJ 


SANDIA  NATIONAL  UBORATORIES 

ATTN.  LI8&SECCLASSIFDIV 

SANDIA  N.ATIONAL  LABORATORIES 
ATTN.  DIV7U1  B  VORTMAN 
ATTN.  J  BANNISTER 
ATTN.  ORG  7111  L  HILL 
ATTN-  TECH  LIB  3141 

OTHER  GOVERNMENT 

FEDERAL  EMERGENCY  MANAGEMENT  AGENCY 

ATTN.  ASST  ASSOC  DIR  FOR  RSCH  J  KERR 
ATTN  W  CHIPMAN/NPCP 

DEPARTMENT  OF  DEFENSE  CONTRACTORS 

AEROSPACE  CORP 

ATTN:  LIBRARY  ACQUISITION  M)/199 

AGBABIAN  ASSOCIATES,  INC 
ATTN-  CBAGGE 

ANALYTIC  SERVICES,  INC  (ANSER) 

ATTN.  G  HESSELBACHER 

APPLIED  RESEARCH  ASSOCIATES,  INC 
2  CYS  ATTN-  B  BINGHAM 
ATTN  N  HIGGINS 

APPLIED  RESEARCH  ASSOCIATES,  INC 
ATTN:  DPIEPENBURG 

APPLIED  RESEARCH  ASSOCIATES,  INC 
ATTN.  P  FRANK 

APPLIED  THEORY,  INC 

ATTN.  JTRUIIG 

AVCO  SYSTEMS  DIVISION 


ATTN 

LIBRARY  A830 

BDM  CORP 

ATTN 

A  LAVAGNINO 

ATTN. 

CORPORATE  LIB 

ATTN. 

T  NEIGHBORS 

BDM  CORP 

ATTN 

F  LEECH 

ATTN 

R  HENSLEY 

BOEING  CO 

ATTN 

AEROSPACE  LIB 

ATTN 

M/S  42/37  R  CARLSON 

BOEING  CO 

ATTN 

MS-85-20  D  CHOATE 

CALIFORNIA  RESEARCH  &  TECHNOLOGY,  INC 
ATTN  KKREYENHAGEN 
ATTN  LIBRARY 

CALIFORNIA  RESEARCH  A  TECHNOLOGY,  INC 
ATTN  F SAUER 
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DEPARTMENT  OF  DEFENSE  CONTRACTORS  (CONTINUED) 

UNIVERSITY  OF  NEW  MEXICO 

CALIFORNIA  RESEARCH  &  TECHNOLOGY,  INC 

NEW  MEXICO  ENGINEERING  RESEARCH  INSTITUTE 
ATTN-  NBAUM 

ATTN  TECH  LIB 

CALSPAN  CORP 

PACIFIC-SIERRA  RESEARCH  CORP 

ATTN;  H  BRODE,  CHAIRMAN  SAGE 

ATTN.  LIBRARY 

CARPENTER  RESEARCH  CORP 

PACIFICA  TECHNOLOGY 

ATTN-  R ALLEN 

ATTN;  HJ  CARPENTER 

ATTN;  RBJORK 

DENVER.  UNIVERSITY  OF 

PATEL  ENTERPRISES,  INC 

ATTN-  SEC  OFCR  FOR  J  WISOTSKI 

ATTN  M  PATEL 

EG&G  VIASH  ANALYTICAL  SVCS  CTR,  INC 

PHYSICAL  RESEARCH.  INC 

ATTN  LIBRARY 

ATTN.  WMENDES 

GENERAL  RESEARCH  CORP 

PHYSICS  INTERNATIONAL  CO 

ATTN  TECHNICAL  INFORMATION  OFFICE 

ATTN-  E  MOORE 

ATTN:  LBEHRMANN 

H-TECH  LABS,  INC 

ATTN:  BHARTENBAUM 

R  &  D  ASSOCIATES 

IIT  RESEARCH  INSTITUTE 

ATTN-  AKUHL 

ATTN:  D  SIMONS 

ATTN-  DOCUMENTS  LIBRARY 

ATTN:  J  LEWIS 

INSTITUTE  FOR  DEFENSE  ANALYSES 

ATTN.  TECH  INFO  CTR 

ATTN.  W  WRIGHT 

ATTN  CUSSIFIEO  LIBRARY 

KAMAN  SCIENCES  CORP 

R&DASSOCIATES 

ATTN:  GGANONG 

ATTN.  LMENTE 

ATTN:  LIBRARY 

RAND  CORP 

ATTN.  RRUETENIK 

ATTN  P  DAVIS 

ATTN  WLEE 

KAMAN  SCIENCES  CORP 

RAND  CORP 

ATTN  B BENNETT 

ATTN  r  SHELTON 

ATTN;  LIBRARY 

S-CUBED 

KAMAN  SCIENCES  CORP 

ATTN.  DGRINE 

ATTN  LIBRARY 

ATTN.  E  CONRAD 

ATTN-  TRINEY 

KAMAN  TEMPO 

SCIENCE  &  ENGRG  ASSOC,  INC 

ATTN  DASIAC 

ATTN.  B  CHAMBERS 

KAMAN  TEMPO 

SCIENCE  APPLICATIONS  INTL  CORP 

ATTN.  DASIAC 

ATTN:  TECH  LIB 

LOCKHEED  MISSILES  &  SPACE  CO.  INC 

SCIENCE  APPLICATIONS  INTL  CORP 

ATTN-  J  WEISNER  DEPT  80-82 

ATTN-  D  MAXWELL 

ATTN:  TECH  INFO  CTR  D/COLL,  D/90-11,  B/106 

MARTIN  MARIETTA  CORP 

SCIENCE  APPLICATIONS  INTL  CORP 

ATTN.  J  COCKAYNE 

ATTN.  GFOTIEO 

ATTN:  MKNASEL 

MCDONNELL  DOUGLAS  CORP 

ATTN.  RSIEVERS 

ATTN  WLAYSON 

ATTN-  RHALPRIN 

MCDONNELL  DOUGLAS  CORP 

SCIENCE  APPLICATIONS  INTL  CORP 

ATTN  GBINNINGER 

ATTN:  M  POTTER  MS/35-18 

MERRITT  CASES,  INC 

SOUTHWEST  RESEARCH  INSTITUTE 

ATTN  A  WENZEL 

ATTN-  J  MERRITT 

ATTN-  W  BAKER 

ATTN  LIBRARY 
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DEPARTMENT  OF  DEFENSE  CONTRACTORS  (CONTINUED) 


TERRA  TEK.  INC 


SRI  INTERNATIONAL 

ATTN:  GABRAHAMSON 

STRUCTURAL  MECHANICS  ASSSOC,  INC 
ATTN.  R  KENNEDY 

TELEDYNE  BROWN  ENGINEERING 
ATTN  D  ORMOND 
ATTN.  F  LEOPARD 
ATTN;  J  RAVENSCRAFT 

TRW  ELECTRONICS  &  DEFENSE  SECTOR 
2CYS  ATTN  NLIPNER 

ATTN:  TECH  INFO  CTR 


ATTN:  LIBRARY 

TRW  ELECTRONICS  &  DEFENSE  SECTOR 
ATTN:  EWONG 
ATTN:  GHULCHER 
ATTN.  PDAI 

WEIDLINGER  ASSOC,  CONSULTING  ENGRG 
ATTN:  TDEEVY 

WEIDLINGER  ASSOC.  CONSULTING  ENGRG 
ATTN:  I  SANDLER 
AnN.  M  BARON 

WEIDLINGER  ASSOC,  CONSULTING  ENGRG 
ATTN:  JISENBERG 
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